/* See float_image_interpolate_C1.h */
/* Last edited on 2009-06-03 23:23:30 by stolfi */ 

#define _GNU_SOURCE
#include <limits.h>
#include <float.h>
#include <assert.h>
#include <math.h>

#include <affirm.h>
#include <bool.h>
#include <indexing.h>
#include <float_image.h>

#include <float_image_interpolate_C0.h>

/* INTERNAL PROTOTYPES */

void float_image_interpolate_C0_get_indices_and_weights(double z, int N, ix_reduction_t red, int i[], double w[]);
  /* Computes the indices {i[0..1]} of the pixels needed to
    interpolate a row (or column) of samples at the fractional
    coordinate {z}, and their respective weights {w[0..1]}. Indices {i}
    outside the range {0..N-1} are reduced according to {ix_reduce(i, N,
    red)}. The reduced index may be -1 to denote `no such pixel'. */

void float_image_interpolate_C0_get_samples_and_weights
  ( float_image_t *A,
    int c,
    double x, 
    double y,
    ix_reduction_t red,
    float *p[],
    double wx[],
    double wy[]
  );
  /* Computes the addresses {p[0..3]} of the 2×2 samples
    of channel {c} that are involved in the interpolation of image {A}
    at the point {(x,y)}, and the corresponding weights. The 
    pixel addresses are stored by rows.

    Indices outside the valid ranges are reduced according to
    {ix_reduce(i, N, red)}. If the reduced index is invalid, 
    the corresponding pointer is set to NULL. */

/* IMPLEMENTATIONS */
  
void float_image_interpolate_C0_get_indices_and_weights
  ( double z, 
    int N,
    ix_reduction_t red, 
    int *i0p, 
    int *i1p, 
    double *w0p,
    double *w1p
  )
  {
    /* Adjust {z} so that the integer part is the first pixel: */
    z = z - 0.5;
    
    /* Get the raw index of the first pixel: */
    int iz = (int)floor(z);
    
    /* Get the fraction {f} fromdata pixel 0: */
    double f = z - iz;

    /* Make sure that the fraction is in {[0 _ 1]}: */
    double fz = z - iz;
    if (fz < 0.0) { iz--; fz += 1.0; }
    if (fz > 1.0) { iz++; fz -= 1.0; }

    /* Compute the reduced indices {i0,i1}: */
    int i0 = ix_reduce(iz + 0, N, red);
    int i1 = ix_reduce(iz + 1, N, red);

    /* Compute the interpolation weights {w0,w1}: */
    /* Linear C0 interpolant, unit sum: */
    double gz = 1.0 - fz;
    double w0 = gz;
    double w1 = fz;

    (*i0p) = i0;
    (*i1p) = i1;
    (*w0p) = w0;
    (*w1p) = w1;
  }

void float_image_interpolate_C0_get_data
  ( float_image_t *A,
    int c,
    double x, 
    double y,
    ix_reduction_t red,
    float **p00p,
    float **p10p,
    float **p01p,
    float **p11p,
    double *wx0p,
    double *wx1p,
    double *wy0p, 
    double *wy1p
  )
  {
    !! revise
    int ix0, ix1;       /* Column indices to interpolate. */
    float_image_interpolate_C0_get_indices_and_weights(x, A->sz[1], red, &ix0, &ix1, wx0p, wx1p);
    int dx = ix1 - ix0; /* Usually 1, may be 0. */
    
    int iy0, iy1;       /* Row indices to interpolate. */
    float_image_interpolate_C0_get_indices_and_weights(y, A->sz[2], red, &iy0, &iy1, wy0p, wy1p);
    int dy = iy1 - iy0; /* Usually 1, may be 0. */
    
    /* Get sample addresses: */
    (*p00p) = float_image_get_sample_address(A, c, ix0, iy0);
    (*p10p) = (*p00p) + dx*A->st[1];
    (*p01p) = (*p00p) + dy*A->st[2];
    (*p11p) = (*p10p) + dy*A->st[2];
  }

float float_image_interpolate_C0_sample(float_image_t *A, int c, double x, double y, ix_reduction_t red)
  { !! revise
    /* Get pixel addresses and weights: */
    double wx0, wx1;    /* Column weights. */
    double wy0, wy1;    /* Sample weights. */
    float *p00, *p10, *p01, *p11; /* Sample addresses. */
    float_image_interpolate_C0_get_data
      (A, c, x, y, red, &p00, &p10, &p01, &p11, &wx0, &wx1, &wy0, &wy1);
    
    /* Apply interpolation formula: */
    return ((*p00)*wx0 + (*p10)*wx1)*wy0 + ((*p01)*wx0 + (*p11)*wx1)*wy1;
  }
  
void float_image_interpolate_C0_pixel(float_image_t *A, double x, double y, ix_reduction_t red, float v[])
  { !! revise
    /* Get pixel addresses and weights: */
    double wx0, wx1;    /* Column weights. */
    double wy0, wy1;    /* Sample weights. */
    float *p00, *p10, *p01, *p11; /* Sample addresses. */
    float_image_interpolate_C0_get_data
      (A, 0, x, y, red, &p00, &p10, &p01, &p11, &wx0, &wx1, &wy0, &wy1);
    
    /* Apply interpolation formula to all channels: */
    int NC = A->sz[0];         /* Number of channels. */
    ix_step_t cst = A->st[0];  /* Position increment between channels. */
    int c;
    for (c = 0; c < NC; c++)
      { /* Interpolate: */
        v[c] = ((*p00)*wx0 + (*p10)*wx1)*wy0 + ((*p01)*wx0 + (*p11)*wx1)*wy1;
        /* Advance sample pointers to next channel: */
        p00 += cst; p10 += cst;  p01 += cst;  p11 += cst;
      }
  }