/*
 * This file is part of the ESO SINFONI Pipeline
 * Copyright (C) 2004,2005 European Southern Observatory
 *
 * This program 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 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
/*----------------------------------------------------------------------------

   File name     :    function_1d.c
   Author         :    Nicolas Devillard
   Created on    :    Tue, Sept 23 1997    
   Description    :    1d signal processing related routines    

 ---------------------------------------------------------------------------*/
/*
 $Id: sinfo_function_1d.c,v 1.7 2012-03-02 08:42:20 amodigli Exp $
 $Author: amodigli $
 $Date: 2012-03-02 08:42:20 $
 $Revision: 1.7 $
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
/*----------------------------------------------------------------------------
                                   Includes
 ---------------------------------------------------------------------------*/
#include <math.h>
#include <string.h>
#include "sinfo_function_1d.h"
#include "sinfo_fit_curve.h"
#include "sinfo_median.h"

/**@{*/
/**
 * @defgroup sinfo_function_1d 1d functions
 *
 * TBD
 */


/*----------------------------------------------------------------------------
                                   Defines
 ---------------------------------------------------------------------------*/
/*
 * This parameter sets up the half size of the domain around which a
 * centroid position will be computed.
 */
#define HALF_CENTROID_DOMAIN    5
/*----------------------------------------------------------------------------
                        Private function prototypes
 ---------------------------------------------------------------------------*/
static double * function1d_generate_smooth_kernel(int filt_type, int hw);

static int
function1d_search_value(
                pixelvalue  *   x,
                int             len,
                pixelvalue      key,
                int         *   foundPtr
) ;

/*----------------------------------------------------------------------------
                              Function codes
 ---------------------------------------------------------------------------*/
/**
  @name        sinfo_function1d_new
  @memo        Allocates a new array of pixelvalues.
  @param    nsamples    Number of values to store in the array.
  @return    Pointer to newly allocated array of pixelvalues.
  @doc

  The returned array must be freed using sinfo_function1d_del(), 
  not cpl_free().
  This is in case some further housekeeping attributes are allocated
  together with the object in the future.

  Returns NULL in case of error.

 */

pixelvalue *
sinfo_function1d_new(int nsamples)
{
    if (nsamples<1) return NULL ;
    return cpl_calloc(nsamples, sizeof(pixelvalue)) ;
}


/**
  @name        sinfo_function1d_del
  @memo        Deallocate an array of pixelvalues.
  @param    s    Array to deallocate.
  @return    void
  @doc

  Deallocates an array allocated by sinfo_function1d_new().
 */
void 
sinfo_function1d_del(pixelvalue * s)
{
    if (s)
        cpl_free(s);
    return ;
}

/**
  @name        sinfo_function1d_dup
  @memo        Copy an array of pixelvalues to a new array.
  @param    arr        Array to copy.
  @param    ns        Number of samples in the array.
  @return    Pointer to newly allocated array of pixelvalues.
  @doc

  Creates a new array using sinfo_function1d_new(), with the same number of
  samples as the input signal, then copies over all values from source
  to destination array using memcpy().

  The returned array must be freed using sinfo_function1d_del(), 
  not cpl_free().

 */
pixelvalue * 
sinfo_function1d_dup(pixelvalue * arr, int ns)
{
    pixelvalue    *    n_arr ;

    n_arr = sinfo_function1d_new(ns);
    memcpy(n_arr, arr, ns * sizeof(pixelvalue));
    return n_arr ;
}

/**
  @name        sinfo_function1d_find_centroid
  @memo        Find out a line centroid to subpixel precision.
  @param    line    Array of pixels.
  @param    npix    Number of pixels in the array.
  @return    Centroid position as a double.
  @doc

  The input signal is assumed to be flat almost everywhere, with a
  single peak somewhere around the middle. Other kinds of signals are
  not handled correctly.

  The position of the peak is located to subpixel precision by
  simply weighting positions with pixelvalues.

  TODO: not used
 */
double
sinfo_function1d_find_centroid(
                pixelvalue    *    line,    /*    the input line    */
                int    npix    /*     number of pixels in this line    */
)
{
    pixelvalue    max ;
    double        centroid ;
    double        weights ;
    int            i, maxpos ;

    /*
     * Search for the maximum pixel value on the line
     */

    max = line[0] ;
    maxpos = 0 ;
    for (i=1 ; i<npix ; i++) {
        if (line[i]>max) {
            max = line[i] ;
            maxpos = i ;
        }
    }

    /*
     * The centroid position is the weighted average over the maximum 
     * pixel neighborhood.
     */

    centroid = 0.0 ;
    weights  = 0.0 ;
    for (i=maxpos-HALF_CENTROID_DOMAIN; 
                    i<=maxpos+HALF_CENTROID_DOMAIN; i++) {
        centroid += (double)line[i] * (double)i ;
        weights  += (double)line[i] ;
    }

    centroid /= weights ;
    return centroid ;    
}

/**
  @name        sinfo_function1d_find_locmax
  @memo        Find out a local maximum in a 1d signal around a position.
  @param    line    Array of pixels.
  @param    npix    Number of pixels in the array.
  @param    where    Where to look around.
  @param    hs        Half-size of the search domain.
  @return    Local maximum position as a double.
  @doc

  The closest local maximum to the given position is located to subpixel
  precision. This precision is achieved by simply weighting positions with
  pixelvalues.

  The 'where' parameter indicates where to look for a maximum as an index
  in the array, i.e. it must lie between 0 and npix-1 (inclusive). The 'hs'
  parameter indicates the half-size of the search domain, i.e. if hs=5 a
  local maximum will be searched +/-5 pixels around the requested position.

  Returns a negative value if an error occurred.

  TODO: not used
 */

double
sinfo_function1d_find_locmax(
                pixelvalue    *    line,
                int                npix,
                int                where,
                int                hs
)
{
    pixelvalue    max ;
    double        centroid ;
    double        weights ;
    int            i, maxpos ;


    if ((where<hs) || (where>(npix-hs-1))) {
        return (double)-1.0 ;
    }

    /*
     * Search for the closest local maximal around the requested range.
     */
    max = line[where] ;
    maxpos = where ;
    for (i=-hs ; i<=hs ; i++) {
        if (line[where+i]>max) {
            max = line[where+i] ;
            maxpos = where+i ;
        }
    }

    /*
     * The centroid position is the weighted average over the maximum 
     * pixel neighborhood.
     */

    centroid = 0.0 ;
    weights  = 0.0 ;
    for (i=maxpos-hs; i<=maxpos+hs; i++) {
        centroid += (double)line[i] * (double)i ;
        weights  += (double)line[i] ;
    }
    if (fabs(weights)>1e-6) {
        centroid /= weights ;
    } else {
        centroid = -1.0 ;
    }
    return centroid ;    
}

/**
  @name        sinfo_function1d_filter_lowpass
  @memo        Apply a low-pass filter to a 1d signal.
  @param    input_sig    Input signal
  @param    samples        Number of samples in the signal
  @param    filter_type    Type of filter to use.
  @param    hw            Filter half-width.
  @return    Pointer to newly allocated array of pixels.
  @doc

  This kind of low-pass filtering consists in a convolution with a
  given kernel. The chosen filter type determines the kind of kernel
  to apply for convolution. Possible kernels and associated symbols
  can be found in function_1d.h.

  Smoothing the signal is done by applying this kind of low-pass
  filter several times.

  The returned smooth signal has been allocated using
  sinfo_function1d_new(), it must be freed using sinfo_function1d_del(). The
  returned signal has exactly as many samples as the input signal.

 */

pixelvalue *
sinfo_function1d_filter_lowpass(
                pixelvalue    *    input_sig,
                int                samples,
                int                filter_type,
                int                hw
)
{
    pixelvalue    *    out_sig ;
    int                i, j ;
    double            replace ;
    double        *    kernel ;

    /* allocate output signal */
    out_sig = sinfo_function1d_new(samples);

    /* generate low-pass filter kernel */
    kernel = function1d_generate_smooth_kernel(filter_type, hw) ;

    /* compute sinfo_edge effects for the first hw elements */
    for (i=0 ; i<hw ; i++) {
        replace = 0.0 ;
        for (j=-hw ; j<=hw ; j++) {
            if (i+j<0) {
                replace += kernel[hw+j] * (double)input_sig[0] ;
            } else {
                replace += kernel[hw+j] * (double)input_sig[i+j] ;
            }
        }
        out_sig[i] = (pixelvalue)replace ;
    }

    /* compute sinfo_edge effects for the last hw elements */
    for (i=samples-hw ; i<samples ; i++) {
        replace = 0.0 ;
        for (j=-hw ; j<=hw ; j++) {
            if (i+j>samples-1) {
                replace += kernel[hw+j] * (double)input_sig[samples-1] ;
            } else {
                replace += kernel[hw+j] * (double)input_sig[i+j] ;
            }
        }
        out_sig[i] = (pixelvalue)replace ;
    }

    /* compute all other elements */
    for (i=hw ; i<samples-hw ; i++) {
        replace = 0.0 ;
        for (j=-hw ; j<=hw ; j++) {
            replace += kernel[hw+j] * (double)input_sig[i+j] ;
        }
        out_sig[i] = (pixelvalue)replace ;
    }

    cpl_free(kernel) ;
    return out_sig ;

}

/**
  @name        function1d_generate_smooth_kernel
  @memo        Generate a kernel for smoothing filters (low-pass).
  @param    filt_type    Type of kernel to generate.
  @param    hw            Kernel half-width.
  @return    Pointer to newly allocated kernel.
  @doc

  Supported kernels can be found in function_1d.h

  The returned array has been allocated with cpl_calloc(), it must be
  freed using cpl_free(). The returned array contains 2h+1 doubles, where
  h is the kernel half-width.

 */

static double * 
function1d_generate_smooth_kernel(int filt_type, int hw)
{
    double  *   kernel ;
    double      norm ;
    int         i ;

    kernel = (double*)cpl_calloc(2*hw+1, sizeof(double)) ;

    switch(filt_type) {

    case LOW_PASS_LINEAR:
        for (i=-hw ; i<=hw ; i++) {
            /* flat kernel */
            kernel[hw+i] = 1.0 / (double)(2*hw+1) ;
        }
        break ;

    case LOW_PASS_GAUSSIAN:
        norm = 0.00 ;
        for (i=-hw ; i<=hw ; i++) {
            /* sinfo_gaussian kernel */
            kernel[hw+i] = exp(-(double)(i*i)) ;
            norm += kernel[hw+i] ;
        }
        for (i=0 ; i<2*hw+1 ; i++) {
            kernel[i] /= norm ;
        }
        break ;

    default:
        sinfo_msg_error("unrecognized low pass filter: "
                        "cannot generate kernel") ;
        cpl_free( kernel );
        return (double*)NULL ;
        break ;
    }

    return kernel ;
}

/**
  @name        sinfo_function1d_median_smooth
  @memo        Apply a 1d sinfo_median filter of given half-width.
  @param    list        List of input pixelvalues.
  @param    np            Number of points in the list.
  @param    hw            Filter half-width.
  @return    Pointer to newly allocated array of pixelvalues.
  @doc

  This function applies a sinfo_median smoothing to a given signal and
  returns a newly allocated signal containing a sinfo_median-smoothed
  version of the input. The returned array has exactly as many samples
  as the input array. It has been allocated using sinfo_function1d_new() and
  must be deallocated using sinfo_function1d_del().

  For half-widths of 1,2,3,4, the filtering is optimized for speed.

  TODO: not used
 */

pixelvalue * 
sinfo_function1d_median_smooth(
                pixelvalue * list,
                int          np,
                int             hw)
{
    int             i,j ;
    pixelvalue    *    row ;
    pixelvalue    *     smoothed ;

    /* simply copy first 3 and last 3 items */
    smoothed = sinfo_function1d_new(np);
    for (i=0 ; i<hw ; i++) {
        smoothed[i] = list[i] ;
    }
    for (i=np-hw ; i<np ; i++) {
        smoothed[i] = list[i] ;
    }

    /* sinfo_median filter on all central items */
    row = sinfo_function1d_new(2*hw+1);
    for (i=hw ; i<np-hw ; i++) {
        for (j=-hw ; j<=hw ; j++) {
            row[j+hw] = list[i+j] ;
        }
        smoothed[i] = sinfo_median_pixelvalue(row, 2*hw+1) ; 
    }
    sinfo_function1d_del(row) ;
    return smoothed ;
}

/**
  @name        sinfo_function1d_remove_lowfreq
  @memo        Subtract low-frequency components from a signal.
  @param    pixel_signal    Input signal.
  @param    ns        Number of samples.
  @return    Pointer to newly allocated array of pixelvalues.
  @doc

  The returned signal is such as: out = in - smooth(in).

  The returned array has been allocated using sinfo_function1d_dup(), it
  must be deallocated using sinfo_function1d_del(). The returned array has
  exactly as many elements as the input array.

TODO: not used
 */

#define LOWFREQ_PASSES        5

pixelvalue * 
sinfo_function1d_remove_lowfreq(
                pixelvalue * pixel_signal,
                int             ns)
{
    pixelvalue    *    sig_in ;
    pixelvalue    *    smooth ;
    int                i ;


    /* Apply severe low-pass filter several times */
    sig_in = sinfo_function1d_dup(pixel_signal, ns);
    for (i=0 ; i<LOWFREQ_PASSES ; i++) {
        smooth = sinfo_function1d_filter_lowpass( sig_in, ns, 
                        LOW_PASS_LINEAR, 5);
        cpl_free(sig_in);
        sig_in = smooth ;
    }

    /* Subtract smoothed signal from input signal */
    for (i=0 ; i<ns ; i++) {
        smooth[i] = pixel_signal[i] - smooth[i];
    }
    return smooth ;
}

#undef LOWFREQ_PASSES

/**
  @name        sinfo_function1d_remove_thermalbg
  @memo        Remove thermal background from a signal.
  @param    pixel_signal    Input signal.
  @param    ns        Number of samples in the input signal.
  @return    Pointer to newly allocated array of pixelvalues.
  @doc

  Many assumptions are made about the input signal. What is expected
  is typically the a collapsed image taken in K band, where the
  thermal background is rising as an exponential of the wavelength.

  This function tries to remove the thermal background contribution by
  first estimating it, then interpolating missing background values,
  and finally subtracting it from the input signal.

  The returned array has been allocated using sinfo_function1d_new(), it
  must be freed using sinfo_function1d_del(). The returned array has exactly
  as many samples as the input array.

  TODO: not used
 */

#define SAMPLE_BORDER    10

pixelvalue * 
sinfo_function1d_remove_thermalbg(
                pixelvalue * pixel_signal,
                int             ns)
{
    pixelvalue    *    smooth ;
    int                i ;

    int                nmin ;
    pixelvalue        lef[2], rig[2];

    pixelvalue    *    x,
    *    y,
    *    spl_x,
    *    spl_y ;
    double            med_y ;
    double            avg2med ;



    /* Detect all local minima */
    nmin = 0 ;
    x = sinfo_function1d_new(ns);
    y = sinfo_function1d_new(ns);

    for (i=SAMPLE_BORDER ; i<(ns-SAMPLE_BORDER) ; i++) {
        lef[0] = pixel_signal[i-2];
        lef[1] = pixel_signal[i-1];
        rig[0] = pixel_signal[i+1];
        rig[1] = pixel_signal[i+2];

        if ( (pixel_signal[i] < lef[0]) &&
                        (pixel_signal[i] < lef[1]) &&
                        (pixel_signal[i] < rig[0]) &&
                        (pixel_signal[i] < rig[1])) {
            x[nmin] = (pixelvalue)i ;
            y[nmin] = pixel_signal[i];
            nmin ++ ;
        }
    }


    /* Interpolate linearly missing values */
    spl_x = sinfo_function1d_new(ns);
    spl_y = sinfo_function1d_new(ns);
    for (i=0 ; i<ns ; i++) {
        spl_x[i] = (pixelvalue)i ;
    }
    sinfo_function1d_interpolate_linear(x, y, nmin, spl_x, spl_y, ns);

    sinfo_function1d_del(x) ;
    sinfo_function1d_del(y) ;
    sinfo_function1d_del(spl_x);

    /* Compute sinfo_median and average distance to the sinfo_median */
    med_y = (double)sinfo_median_pixelvalue(pixel_signal, ns);
    avg2med = 0.0 ;
    for (i=0 ; i<ns ; i++) {
        avg2med += fabs((double)pixel_signal[i] - med_y) ;
    }
    avg2med /= (double)ns ;

    /* Reset all pixels out of sinfo_median + 2 * avg2med to zero. */
    for (i=0 ; i<ns ; i++) {
        double dist = fabs((double)pixel_signal[i] - med_y);
        if (dist > (2.0*avg2med)) {
            spl_y[i] = (pixelvalue)0 ;
        }
    }


    smooth = sinfo_function1d_new(ns);
    for (i=0 ; i<ns ; i++) {
        if (spl_y[i]>1e-4) {
            smooth[i] = pixel_signal[i] - spl_y[i];
        } else {
            smooth[i] = 0.0 ;
        }
    }
    sinfo_function1d_del(spl_y);
    return smooth ;
}

#undef LOWFREQ_PASSES

/**
  @name        sinfo_function1d_interpolate_linear
  @memo        Linear signal interpolation.
  @param    x        Input list of x positions.
  @param    y        Input list of y positions.
  @param    len        Number of samples in x and y.
  @param    spl_x    List of abscissas where the signal must be computed.
  @param    spl_y    Output list of computed signal values.
  @param    spl_len    Number of samples in spl_x and spl_y.
  @return    void
  @doc

  To apply this interpolation, you need to provide a list of x and y
  positions, and a list of x positions where you want y to be computed
  (with linear interpolation).

  The returned signal has spl_len samples. It has been allocated using
  sinfo_function1d_new() and must be deallocated using sinfo_function1d_del().

 */

void 
sinfo_function1d_interpolate_linear(
                pixelvalue    *    x,
                pixelvalue    *    y,
                int                len,
                pixelvalue    *    spl_x,
                pixelvalue    *    spl_y,
                int                spl_len
)
{
    double        a, b ;
    int            i, j ;


    for (i=0 ; i<spl_len ; i++) {
        /* Find (x1,y1) on the left of the current point */
        int found = 0 ;
        for (j=0 ; j<(len-1) ; j++) {
            if ((spl_x[i]>=x[j]) && (spl_x[i]<=x[j+1])) {
                found++ ;
                break ;
            }
        }
        if (!found) {
            spl_y[i] = 0.0;
        } else {
            a = ((double)y[j+1]-(double)y[j]) /
                            ((double)x[j+1]-(double)x[j]);
            b = (double)y[j] - a * (double)x[j] ;
            spl_y[i] = (pixelvalue)(a * (double)spl_x[i] + b) ;
        }
    }
    return ;
}

/**
  @name        function1d_search_value
  @memo        Conducts a binary search for a value.
  @param    x     Contains the abscissas of interpolation.
  @param    len     Length of the x array.
  @param    key     The value to locate in x.
  @param    foundPtr Output flag, 1 if value was found, else 0.
  @return    The index of the largest value in x for which x[i]<key.
  @doc

  This function does a binary search for a value in an array. This
  routine is to be called only if key is in the interval between x[0]
  and x[len-1]. The input x array is supposed sorted.

 */

static int
function1d_search_value(
                pixelvalue    *    x,
                int             len,
                pixelvalue         key,
                int         *    foundPtr
)
{


    int low  = 0;
    int high = len - 1;

    while (high >= low) {
        int middle = (high + low) / 2;
        if (key > x[middle]) {
            low = middle + 1;
        } else if (key < x[middle]) {
            high = middle - 1;
        } else {
            *foundPtr = 1;
            return (middle);
        }
    }
    *foundPtr = 0;
    return (low);
}

/**
  @name        sinfo_function1d_natural_spline
  @memo        Interpolate a sinfo_vector along new abscissas.
  @param    x        List of x positions.
  @param    y        List of y positions.
  @param    len        Number of samples in x and y.
  @param    splX    Input new list of x positions.
  @param    splY    Output list of interpolated y positions.
  @param    splLen    Number of samples in splX and splY.
  @return    Int 0 if Ok, -1 if error.
  @doc

  Reference:

  \begin{verbatim}
      Numerical Analysis, R. Burden, J. Faires and A. Reynolds.
      Prindle, Weber & Schmidt 1981 pp 112
  \end{verbatim}

  Provide in input a known list of x and y values, and a list where
  you want the signal to be interpolated. The returned signal is
  written into splY.

 */

int
sinfo_function1d_natural_spline(
                pixelvalue    *     x,
                pixelvalue    *     y,
                int             len,
                pixelvalue    *     splX,
                pixelvalue    *     splY,
                int             splLen
)
{
    int             end;
    int             loc,
    found;
    register int     i,
    j,
    n;
    double         *    h; /* sinfo_vector of deltas in x */
    double         *    alpha;
    double         *    l,
    *    mu,
    *    z,
    *    a,
    *    b,
    *    c,
    *    d;

    end = len - 1;

    a = cpl_malloc(sizeof(double) * splLen * 9) ;
    b = a + len;
    c = b + len;
    d = c + len;
    h = d + len;
    l = h + len;
    z = l + len;
    mu = z + len;
    alpha = mu + len;

    for (i = 0; i < len; i++) {
        a[i] = (double)y[i];
    }

    /* Calculate sinfo_vector of differences */
    for (i = 0; i < end; i++) {
        h[i] = (double)x[i + 1] - (double)x[i];
        if (h[i] < 0.0) {
            cpl_free(a) ;
            return -1;
        }
    }

    /* Calculate alpha sinfo_vector */
    for (n = 0, i = 1; i < end; i++, n++) {
        /* n = i - 1 */
        alpha[i] = 3.0 * ((a[i+1] / h[i]) - (a[i] / h[n]) - (a[i] / h[i]) +
                        (a[n] / h[n]));
    }

    /* Vectors to solve the tridiagonal sinfo_matrix */
    l[0] = l[end] = 1.0;
    mu[0] = mu[end] = 0.0;
    z[0] = z[end] = 0.0;
    c[0] = c[end] = 0.0;

    /* Calculate the intermediate results */
    for (n = 0, i = 1; i < end; i++, n++) {
        /* n = i-1 */
        l[i] = 2 * (h[i] + h[n]) - h[n] * mu[n];
        mu[i] = h[i] / l[i];
        z[i] = (alpha[i] - h[n] * z[n]) / l[i];
    }
    for (n = end, j = end - 1; j >= 0; j--, n--) {
        /* n = j + 1 */
        c[j] = z[j] - mu[j] * c[n];
        b[j] = (a[n] - a[j]) / h[j] - h[j] * (c[n] + 2.0 * c[j]) / 3.0;
        d[j] = (c[n] - c[j]) / (3.0 * h[j]);
    }

    /* Now calculate the new values */
    for (j = 0; j < splLen; j++) {
        double v = (double)splX[j];
        splY[j] = (pixelvalue)0;

        /* Is it outside the interval? */
        if ((v < (double)x[0]) || (v > (double)x[end])) {
            continue;
        }
        /* Search for the interval containing v in the x sinfo_vector */
        loc = function1d_search_value(x, len, (pixelvalue)v, &found);
        if (found) {
            splY[j] = y[loc];
        } else {
            loc--;
            v -= (double)x[loc];
            splY[j] = (pixelvalue)(    a[loc] +
                            v * (b[loc] +
                                            v * (c[loc] +
                                                            v * d[loc])));
        }
    }
    cpl_free(a) ;
    return 0;
}

/**
  @name        sinfo_function1d_average_reject
  @memo        Sorts the input signal, takes out highest and lowest
          values, and returns the average of the remaining pixels.
  @param    line        Input signal.
  @param    npix        Number of samples in the input signal.
  @param    pix_low        Number of lowest pixels to reject.
  @param    pix_high    Number of highest pixels to reject.
  @return    The filtered average of the input signal.
  @doc

  No input parameter is modified.

  The input signal is first copied. This copy is then sorted, and the
  highest and lowest pixels are taken out of the list. Remaining
  pixelvalues are averaged and the result is returned.

  TODO: not used
 */

pixelvalue
sinfo_function1d_average_reject(
                pixelvalue    *    line,
                int                npix,
                int                pix_low,
                int                pix_high)
{
    pixelvalue    *    sorted ;
    int                i ;
    double            avg ;

    /* Sanity tests */
    if ((line==NULL) || (npix<1)) return (pixelvalue)0 ;
    if ((pix_low+pix_high)>=npix) return (pixelvalue)0 ;

    /* Copy input line and sort it */
    sorted = cpl_malloc(npix * sizeof(pixelvalue)) ;
    memcpy(sorted, line, npix * sizeof(pixelvalue)) ;
    sinfo_pixel_qsort(sorted, npix);

    /* Find out average of remaining values */
    avg = 0.00 ;
    for (i=pix_low+1 ; i<(npix-pix_high) ; i++) {
        avg += (double)sorted[i] ;
    }
    cpl_free(sorted);
    avg /= (double)(npix - pix_high - pix_low) ;

    return (pixelvalue)avg ;
}

/**
  @name        sinfo_function1d_xcorrelate
  @memo        Cross-sinfo_correlation of two 1d signals.
  @param    line_i        The reference signal.
  @param    width_i        Number of samples in reference signal.
  @param    line_t        Candidate signal to compare.
  @param    width_t        Number of samples in candidate signal.
  @param    half_search    Half-size of the search domain.
  @param    delta        Output sinfo_correlation offset.
  @return    Maximum cross-sinfo_correlation value as a double.
  @doc

  Two signals are expected in input of this function: a reference
  signal and a candidate signal. They are expected to be roughly the
  same signal up to an offset.

  A cross-sinfo_correlation is computed on 2*half_search+1 values. The
  maximum of likelihood is the maximum cross-sinfo_correlation value between
  signals. The offset corresponding to this position is returned.

  Returns -100.0 in case of error. Normally, the cross-sinfo_correlation
  coefficient is normalized so it should stay between -1 and +1.

  TODO: not used
 */


#define STEP_MIN    (-half_search)
#define STEP_MAX    (half_search)

double 
sinfo_function1d_xcorrelate(
                pixelvalue *    line_i,
                int             width_i,
                pixelvalue *    line_t,
                int             width_t,
                int                half_search,
                double     *    delta
)
{
    double * xcorr ;
    double   xcorr_max ;
    double   mean_i, mean_t ;
    double   rms_i, rms_t ;
    double   sum, sqsum ;
    double   norm ;
    int      maxpos ;
    int      nsteps ;
    int      i ;
    int      step ;



    /* Compute normalization factors */
    sum = sqsum = 0.00 ;
    for (i=0 ; i<width_i ; i++) {
        sum += (double)line_i[i] ;
        sqsum += (double)line_i[i] * (double)line_i[i];
    }
    mean_i = sum / (double)width_i ;
    sqsum /= (double)width_i ;
    rms_i = sqsum - mean_i*mean_i ;

    sum = sqsum = 0.00 ;
    for (i=0 ; i<width_t ; i++) {
        sum += (double)line_t[i] ;
        sqsum += (double)line_t[i] * (double)line_t[i];
    }
    mean_t = sum / (double)width_t ;
    sqsum /= (double)width_t ;
    rms_t = sqsum - mean_t*mean_t ;

    norm = 1.00 / sqrt(rms_i * rms_t);

    nsteps = (STEP_MAX - STEP_MIN) +1 ;
    xcorr = cpl_malloc(nsteps * sizeof(double));
    for (step=STEP_MIN ; step<=STEP_MAX ; step++) {
        xcorr[step-STEP_MIN] = 0.00 ;
        int nval = 0 ;
        for (i=0 ; i<width_t ; i++) {
            if ((i+step > 0) &&
                            (i+step < width_i)) {
                xcorr[step-STEP_MIN] += ((double)line_t[i] - mean_t) *
                                ((double)line_i[i+step] - mean_i) *
                                norm ;
                nval++ ;
            }
        }
        xcorr[step-STEP_MIN] /= (double)nval ;
    }
    xcorr_max = xcorr[0] ;
    maxpos    = 0 ;
    for (i=0 ; i<nsteps ; i++) {
        if (xcorr[i]>xcorr_max) {
            maxpos = i ;
            xcorr_max = xcorr[i];
        }
    }
    cpl_free(xcorr);
    (*delta) = + ((double)STEP_MIN + (double)maxpos);
    return xcorr_max ;
}
/**@}*/