filtering.cpp

#include<cmath>
#include<iostream>
#include<cassert>
#include<opencv/cv.h>
extern "C" {
#include<luaT.h>
#include<TH/TH.h>
}
using namespace std;

typedef THFloatTensor Tensor;
const char* idreal = "torch.FloatTensor";
const char* idfloat = "torch.FloatTensor";
const char* idlong = "torch.LongTensor";
#define Tensor_(a) THFloatTensor_##a
typedef float Real;
typedef double accreal;
typedef unsigned char byte;
typedef unsigned short uint16;
typedef cv::Mat_<float> matf;

static int IntegralImage(lua_State *L) {
  THFloatTensor* input = (THFloatTensor*)luaT_checkudata(L, 1, idfloat);
  THFloatTensor* output = (THFloatTensor*)luaT_checkudata(L, 2, idfloat);

  const int h = input->size[0];
  const int w = input->size[1];
  float* ip = THFloatTensor_data(input);
  float* op = THFloatTensor_data(output);
  const long* const is = input->stride;
  const long* const os = output->stride;

  const matf input_cv(h, w, ip, is[0]*sizeof(float));
  matf output_cv(h+1, w+1, op, os[0]*sizeof(float));
  cv::integral(input_cv, output_cv, CV_32F);

  return 0;
}

inline float iimageSum(const float* iimage, long is,
		       int x1, int y1, int x2, int y2) {
  return iimage[is*y1+x1]+iimage[is*y2+x2]
    - iimage[is*y1+x2] - iimage[is*y2+x1];
}

static int FilterImage(lua_State *L) {
  THFloatTensor* iimage  = (THFloatTensor*)luaT_checkudata(L, 1, idfloat);
  THLongTensor * filters = (THLongTensor *)luaT_checkudata(L, 2, idlong );
  THLongTensor * output  = (THLongTensor *)luaT_checkudata(L, 3, idlong );
  int hmax = lua_tointeger(L, 4);
  int wmax = lua_tointeger(L, 5);

  const int h = iimage->size[0]-1;
  const int w = iimage->size[1]-1;
  const int N = filters->size[0];
  const float* const ip = THFloatTensor_data(iimage );
  const long*  const fp = THLongTensor_data (filters);
  long* const        op = THLongTensor_data (output );
  const long* const is = iimage->stride;
  const long* const fs = filters->stride;
  const long* const os = output->stride;
  assert(is[1] == 1);

  //printf("h=%d\tw=%d\tN=%d\n",h, w, N);
  //printf("is=%d\tfs=%d\tos=%d\tos1=%d\n",is[0], fs[0], os[0], os[1]);
  //printf("hmax=%d\twmax=%d\n",hmax,wmax);

  int i, j, k;
  const long* fp0;
  long* op0;
#pragma omp parallel for private(i, j, k, fp0, op0)
  for (i = 0; i < h-hmax; ++i) {
    op0 = op + i*os[0];
    for (j = 0; j < w-wmax; ++j, op0 += os[1])
      for (fp0 = fp, k=0; k < N; fp0 += fs[0], ++k) {
	const int bit = iimageSum(ip, is[0], j+fp0[1],i+fp0[0], j+fp0[3],i+fp0[2])
	  > iimageSum(ip, is[0], j+fp0[5],i+fp0[4], j+fp0[7],i+fp0[6]);
        *op0 |= bit << k;
      }
  }

  return 0;
}

static const struct luaL_reg libfiltering[] = {
  {"integralImage", IntegralImage},
  {"filterImage", FilterImage},
  {NULL, NULL}
};

LUA_EXTERNC int luaopen_libfiltering (lua_State *L) {
  luaL_openlib(L, "libfiltering", libfiltering, 0);
  return 1;
}