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, int 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){
    int 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;
    int 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;
}