interpolation-test.cpp

#include <assert.h>

#include <inttypes.h>

#include <iostream>
#include <iterator>
#include <vector>
#include <chrono>


template <typename MapT, typename ValueT>
class Interpolator {
public:
  
  Interpolator( const MapT maps[], const ValueT values[], size_t len) {
    assert(len > 0);
    map_.assign(maps, maps + len);
    values_.assign(values, values + len);
  }
  
  ValueT interpolate (const MapT& x) {

    ValueT val;
    if (x < map_.front()) {
      val = values_.front();
    } else if (x > map_.back()) {
      val = values_.back();
    } else {
      // find bucket where the input value falls
      typename std::vector<MapT>::iterator it = std::lower_bound(map_.begin(), map_.end(), x);

      // figure out index into map array to use in values array
      size_t offset = std::distance( map_.begin(), it);
      
      if (it == map_.begin()) {
	val =  values_.front();
      } else {
	
	// now get corresponding output values
	MapT x1 = *it;
	--it;
	MapT x0 = *it;
	
	// now get corresponding output values
	typename std::vector<ValueT>::iterator vit = values_.begin() + offset;
	ValueT y1 = *vit;
	--vit;
	ValueT y0 = *vit;
	ValueT p = static_cast<ValueT>(x - x0) / static_cast<ValueT>(x1 - x0);
	val = (1.0 - p) * y0 + p * y1;
      }
    }
    return val;
  }
  
private:
  std::vector<MapT> map_;
  std::vector<ValueT> values_;
  
};

static const size_t MAP_SIZE = 30;     // size zoom to magnification array
  

static const uint16_t raw_zoom[MAP_SIZE] = {0x0000,
					    0x2372,
					    0x3291,
					    0x3b83,
					    0x41B0,
					    0x4668,
					    0x49FB,
					    0x4D3C,
					    0x5000,
					    0x5270,
					    0x548D,
					    0x56AA,
					    0x589E,
					    0x5A68,
					    0x5BD3,
					    0x5D2B,
					    0x5E4F,
					    0x5F48,
					    0x6018,
					    0x60BF,
					    0x6165,
					    0x61E2,
					    0x625F,
					    0x62B2,
					    0x6306,
					    0x6359,
					    0x6383,
					    0x63AC,
					    0x63D6,
					    0x6400 };
  
  
static  const double mag_out[MAP_SIZE] = {1.0,
					 2.0,
					 3.0,
					 4.0,
					 5.0,
					 6.0,
					 7.0,
					 8.0,
					 9.0,
					 10.0,
					 11.0,
					 12.0,
					 13.0,
					 14.0,
					 15.0,
					 16.0,
					 17.0,
					 18.0,
					 19.0,
					 20.0,
					 21.0,
					 22.0,
					 23.0,
					 24.0,
					 25.0,
					 26.0,
					 27.0,
					 28.0,
					 29.0,
					 30.0 };
  
float getZoom(uint16_t zoom_in)
{
  
  if (zoom_in <= raw_zoom[0]) {
    return mag_out[0];
  }
  
  if (zoom_in > raw_zoom[MAP_SIZE - 1]) {
    return mag_out[MAP_SIZE - 1];
  }
  
  int ii = 0;
  for(; ii < MAP_SIZE - 1; ii++) {
    if ((zoom_in >= raw_zoom[ii]) && (zoom_in < raw_zoom[ii + 1])) {
      break;
    }
  }
  
  double x0 = raw_zoom[ii];
  double x1 = raw_zoom[ii + 1];
  double y0 = mag_out[ii];
  double y1 = mag_out[ii + 1];
  
  double a = (y1 - y0) / (x1 - x0);
  double b = -a * x0 + y0;
  double y = a * zoom_in + b;
  return static_cast<float>(y);
}

void test ()
{
  Interpolator<uint16_t, double> interp( raw_zoom, mag_out, MAP_SIZE);
  uint16_t vals[] = {0x0000,
		     0x1234, 0x2345, 0x3456, 0x4567,
		     0x4FFF, 0x5000, 0x5001,
		     0x589E,
		     0x63FF, 0x6400, 0x6401};
  size_t nvals = sizeof(vals)/sizeof(*vals);
  std::cout << "sizeof(vals) is " << nvals << std::endl;
  for(int ii=0; ii < nvals; ++ii) {
    uint16_t v = vals[ii];
    float zoom = getZoom( v );
    float interp_zoom = static_cast<float>(interp.interpolate( v ));
    std::cout << "for raw mag value 0x" << std::hex << v << " zoom=" << zoom 
	      << " and interp_zoom=" << interp_zoom << std::endl;
  }

  std::chrono::high_resolution_clock::time_point t1, t2;
  double duration;
  
  t1 = std::chrono::high_resolution_clock::now();
  float total_zoom = 0;
  for(int jj = 0; jj < 0x7000; jj++) {
    float zoom = getZoom( jj );
    total_zoom += zoom;
  }
  t2 = std::chrono::high_resolution_clock::now();
  duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
  std::cout << "zoom runtime in microseconds:  " << duration << std::endl;
  std::cout << "total zoom " << total_zoom << std::endl;
  
  t1 = std::chrono::high_resolution_clock::now();
  float total_interp_zoom = 0;
  for(int jj = 0; jj < 0x7000; jj++) {
    float zoom = interp.interpolate( jj );
    total_interp_zoom += zoom;
  }
  t2 = std::chrono::high_resolution_clock::now();
  duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
  std::cout << "interp zoom runtime in microseconds:  " << duration << std::endl;
  std::cout << "total interp zoom " << total_interp_zoom << std::endl;


}


int main(int argc, char** argv) {
  test();
  return 0;
}

run-results.sh

g++ -o interp interp.cpp -lstdc++
Stoll-MBP:cyphy pstoll$ ./interp
sizeof(vals) is 12
for raw mag value 0x0 zoom=1 and interp_zoom=1
for raw mag value 0x1234 zoom=1.51356 and interp_zoom=1.51356
for raw mag value 0x2345 zoom=1.99504 and interp_zoom=1.99504
for raw mag value 0x3456 zoom=3.19782 and interp_zoom=3.19782
for raw mag value 0x4567 zoom=5.78725 and interp_zoom=5.78725
for raw mag value 0x4fff zoom=8.99859 and interp_zoom=8.99859
for raw mag value 0x5000 zoom=9 and interp_zoom=9
for raw mag value 0x5001 zoom=9.0016 and interp_zoom=9.0016
for raw mag value 0x589e zoom=13 and interp_zoom=13
for raw mag value 0x63ff zoom=29.9762 and interp_zoom=29.9762
for raw mag value 0x6400 zoom=30 and interp_zoom=30
for raw mag value 0x6401 zoom=30 and interp_zoom=30
zoom runtime in microseconds:  2345
total zoom 232264
interp zoom runtime in microseconds:  13978
total interp zoom 232264