/************************
* This routine provides interpolation for IR images
* Interpolation is done as nearest neighbour in the vertical
* direction only. In other words, we only interpolate in the
* spacial direction.
* This is written in C for efficiency issues and should
* be loaded as an IDL DLM
*
* IDL interface:
* result = DICal_IRInterp(imgIn, toInterp, badPixs)
*
* INPUTS:
* imgIn - Image array to work on
* toInterp - An array of indicies of pixels to be interpolated over. Must be
*            ordered
* badPixs - An array of indicies of invalid pixels. Must be ordered
*
* RETURNS:
* A double array with the pixels interpolated over
*
* MODIFICATIONS:
*  Jan 30, 2005 Mark Desnoyer - Created
*
*************************/
#include "dical_irinterp.h"
#include "dical_interptools.h"

// Actual Interpolation routine
IDL_VPTR DICal_IRInterp(int argc, IDL_VPTR argv[])
{
  IDL_VPTR imgIn, imgOut, toInterp, badPixs;
  double *in, *out;
  IDL_LONG *toDo, *bad;
  int imgWid, doCnt, badCnt;
  int i,j,n,m;
  int row;
  int found;

  // Allocate the input variables
  imgIn = imgOut = argv[0];
  toInterp = argv[1];
  badPixs = argv[2];

  if(interpCheckParams(imgIn, &imgOut, toInterp, badPixs, &in, &out, &imgWid,
  &toDo, &doCnt, &bad, &badCnt))
    return imgIn;

  // Create the output image
  m = imgOut->value.arr->n_elts;
  for (i=0;i<m;i++){
    // If we need to interpolate
    if (*toDo == i){
      toDo++;
      out[i] = 0;
      // Search vertically from the pixel for the nearest valid neighbor
      row = i/imgWid;
      n = ((row > (imgIn->value.arr->n_elts/imgWid)-row)?
        row:(imgIn->value.arr->n_elts/imgWid)-row)-1;
      found = 0;
      for (j=1;j<n && !found;j++){
        if (isValidPix(bad,badCnt,i+imgWid*j, m-1)){
          out[i] = in[i+imgWid*j];
          found = 1;
        }
        if (isValidPix(bad,badCnt,i-imgWid*j, m-1)){
          out[i] = (out[i]+in[i-imgWid*j])/(++found);
        }
      }

    } else // Copy the value
      out[i] = in[i];
  }

  return imgOut;
}