sicspsi/rebin.c

/**
 * This is a couple of functions to do rebinning. This is neeeded when 
 * you have a number of data points with random coordinates and you wish to 
 * get this back onto a regular grid for better handling. 
 * 
 * copyright: see file COPYRIGHT
 * 
 * Mark Koennecke, May 2007
 */
#include <math.h>
#include "rebin.h"
/*---------------------------------------------------------------------------*/
static int checkPoint(pNXDS target, double x, double y)
{
  int maxx, maxy;

  maxx = getNXDatasetDim(target, 0);
  maxy = getNXDatasetDim(target, 1);

  if (x < .0 || y < .0 || x > maxx - 1 || y > maxy - 1) {
    return 0;
  }
  return 1;
}

/*--------------------------------------------------------------------------*/
static void addNXDatasetValue(pNXDS target, int64_t iPos[], double value)
{
  double val;

  val = value + getNXDatasetValue(target, iPos);
  putNXDatasetValue(target, iPos, val);
}

/*---------------------------------------------------------------------------*/
static void roundPoint(pNXDS target, double x, double y, double value)
{
  int64_t iDim[2];

  iDim[0] = (int) floor(x + .5);
  iDim[1] = (int) floor(y + 0.5);

  addNXDatasetValue(target, iDim, value);
}

/*----------------------------------------------------------------------------
 * I am using the squares here, no need to do the sqrt which would give the 
 * actual distance.
 * ---------------------------------------------------------------------------*/
static double dist(double x1, double y1, double x2, double y2)
{
  double dx, dy;

  dx = x1 - x2;
  dy = y1 - y2;
  return dx * dx + dy * dy;
}

/*-----------------------------------------------------------------------------
 * This tries to distribute value over all of the four neighbouring points 
 * weighted with the distance to them.
-------------------------------------------------------------------------------*/
static void distribute(pNXDS target, double x, double y, double value)
{
  double lldist, lrdist, uldist, urdist, totalDist, ix, iy, frac;
  int64_t iPos[2];

  ix = floor(x);
  iy = floor(y);
  lldist = dist(x, y, ix, iy);
  if (lldist < .01) {
    iPos[0] = ix;
    iPos[1] = iy;
    addNXDatasetValue(target, iPos, value);
    return;
  }

  ix += 1;
  lrdist = dist(x, y, ix, iy);
  if (lrdist < .01) {
    iPos[0] = ix;
    iPos[1] = iy;
    addNXDatasetValue(target, iPos, value);
    return;
  }

  iy += 1;
  urdist = dist(x, y, ix, iy);
  if (urdist < .01) {
    iPos[0] = ix;
    iPos[1] = iy;
    addNXDatasetValue(target, iPos, value);
    return;
  }

  ix = floor(x);
  uldist = dist(x, y, ix, iy);
  if (uldist < .01) {
    iPos[0] = ix;
    iPos[1] = iy;
    addNXDatasetValue(target, iPos, value);
    return;
  }

  lldist = 1. / lldist;
  lrdist = 1. / lrdist;
  urdist = 1. / urdist;
  uldist = 1. / uldist;
  totalDist = lldist + lrdist + urdist + uldist;

  if (totalDist < .01) {
    iPos[0] = (int) floor(x);
    iPos[1] = (int) floor(y);
    addNXDatasetValue(target, iPos, value);
    return;
  }

  iPos[0] = (int) floor(x);
  iPos[1] = (int) floor(y);
  frac = (lldist / totalDist) * value;
  addNXDatasetValue(target, iPos, frac);

  iPos[0] += 1;
  frac = (lrdist / totalDist) * value;
  addNXDatasetValue(target, iPos, frac);

  iPos[1] += 1;
  frac = (urdist / totalDist) * value;
  addNXDatasetValue(target, iPos, frac);

  iPos[0] = (int) floor(x);
  frac = (uldist / totalDist) * value;
  addNXDatasetValue(target, iPos, frac);
}

/*---------------------------------------------------------------------------*/
int rebinPoint2D(pNXDS target, double x, double y, double value,
                 int method)
{
  if (!checkPoint(target, x, y)) {
    return 0;
  }
  switch (method) {
  case TRILINEAR:
    distribute(target, x, y, value);
    break;
  case ROUND:
    roundPoint(target, x, y, value);
    break;
  default:
    return 0;
  }
  return 1;
}