// BLInterp.h
//	Copyright 1994, Center for Computer Vision and Visualization,
//	University of Florida.  All rights reserved.
//

#ifndef BLInterp_h_
#define BLInterp_h_

#include	<math.h>

#include "UcharDI.h"
#include "IntDI.h"
#include "FloatDI.h"
#include "FloatCI.h"
#include "ClosureI.h"

template <class T>
class IA_BLInterp : public IA_ClosureI<IA_Point<double>,T> {
    IA_Image<IA_Point<int>,T>	im;
public:
    IA_BLInterp(IA_Image<IA_Point<int>,T> im_)
    :IA_ClosureI<IA_Point<double>,T>(
	IA_universal_fpset(im_.domain().dim())),
     im(im_) {}

    T operator()(const IA_Point<double> &fp) const;
    IA_ClosureI<IA_Point<double>,T> *clone_self() const {
	return new IA_BLInterp(*this); }
};


template <class T>
T IA_BLInterp<T>::operator() (const IA_Point<double> &fp) const
{
    int      ncorners = 1<<fp.dim();
    T       *vals     = new T[ncorners];
    int     *point    = new int[fp.dim()];
    double  *offset   = new double[fp.dim()];

    int i;
    for (i=0; i<fp.dim(); i++)
	offset[i] = fp[i] - floor(fp[i]);

    for (i=0; i<ncorners; i++) {
	int	k=i;
	for (int idx=0; idx<fp.dim(); idx++) {
	    point[idx] = (int)floor(fp[idx]) + (k&1);
	    k >>= 1;
	}
	IA_Point<int>  p(point, fp.dim());
	vals[i] = im.domain().contains(p) ? im(p) : 0;
    }
    
    int	k=fp.dim()-1;
    ncorners >>=1;
    while (ncorners) {
	for (i=0; i<ncorners; i++) {
	    vals[i] = (1-offset[k])*vals[i] + offset[k]*vals[i+ncorners];
	}
	ncorners >>=1;
	k--;
    }

    T	rval = vals[0];
    delete[] offset;
    delete[] vals;
    delete[] point;

    return rval;
}

#endif
