oclimgutil.cl

// Copyright Naoki Shibata 2018. Distributed under the MIT License.

#ifndef NDEBUG
#define ASSERT(x, str) { if (!(x)) { printf((constant char *) "Assertion failure : %s\n", str); } }
#else
#define ASSERT(x, str) {}
#endif

typedef char int8_t;
typedef short int16_t;
typedef int int32_t;
typedef long int64_t;
typedef unsigned char uint8_t;
typedef unsigned short uint16_t;
typedef unsigned int uint32_t;
typedef unsigned long uint64_t;

//

extern __constant unsigned short s2l[];
extern __constant uchar l2s[];
extern __constant unsigned short cfunc[];
extern __constant unsigned short cfunc2[];
extern __constant float iircoef[32][15];

//

uint32_t packlab(float3 clab) {
  uint32_t ret = 0;
  ret = clamp(convert_uint_rtn(clab.z * 1024), 0U, 1023U);
  ret = (ret << 10) | clamp(convert_uint_rtn(clab.y * 1024), 0U, 1023U);
  ret = (ret << 12) | clamp(convert_uint_rtn(clab.x * 4096), 0U, 4095U);
  return ret;
}

float3 unpacklab(uint32_t plab) {
  float3 f = convert_float3((int3) (plab, plab >> 12, plab >> 22) & (int3) (4095, 1023, 1023));
  return f * (float3) (1.0f/4096, 1.0f/1024, 1.0f/1024) + (float3)(0.5f/4096, 0.5f/1024, 0.5f/1024);
}

int mirror(int x, int y, int iw, int ih) {
  int2 xy = (int2)(x, y);
  xy = clamp(xy, -xy, (int2)(iw, ih) * 2 - 2 - xy);
  return xy.x + xy.y * iw;
}

int mirror1(int x, int iw) {
  return clamp(x, -x, iw * 2 - 2 - x);
}

int repeat1(int x, int iw) {
  x = x < 0 ? x+iw : x;
  x = x >= iw ? x-iw : x;
  return x;
}

int repeat(int x, int y, int iw, int ih) {
  x = x < 0 ? x+iw : x;
  x = x >= iw ? x-iw : x;
  y = y < 0 ? y+ih : y;
  y = y >= ih ? y-ih : y;
  return iw * y + x;
}

float bicubicSub(float p0, float p1, float p2, float p3, float x) {
  float u, v, w;
  v = p1 - p2;
  w = p3 - p0;
  u = v * 3.0f + w;
  u = u * x   + (-4.0f * v + (p0 - p1 - w));
  u = u * x   + (p2 - p0);
  u = u * x * 0.5f + p1;
  return u;
}

float3 bicubicSub3(float3 p0, float3 p1, float3 p2, float3 p3, float x) {
  float3 u, v, w;
  v = p1 - p2;
  w = p3 - p0;
  u = v * 3.0f + w;
  u = u * x    + (-4.0f * v + (p0 - p1 - w));
  u = u * x    + (p2 - p0);
  u = u * x * 0.5f + p1;
  return u;
}

float bicubic(global float *p, float x, float y, int iw, int ih) {
  const int ix = (int)x, iy = (int)y;

  return bicubicSub(bicubicSub(p[mirror(ix-1, iy-1,iw,ih)], p[mirror(ix  , iy-1,iw,ih)], p[mirror(ix+1, iy-1,iw,ih)], p[mirror(ix+2, iy-1,iw,ih)], x-ix),
		    bicubicSub(p[mirror(ix-1, iy  ,iw,ih)], p[mirror(ix  , iy  ,iw,ih)], p[mirror(ix+1, iy  ,iw,ih)], p[mirror(ix+2, iy  ,iw,ih)], x-ix),
		    bicubicSub(p[mirror(ix-1, iy+1,iw,ih)], p[mirror(ix  , iy+1,iw,ih)], p[mirror(ix+1, iy+1,iw,ih)], p[mirror(ix+2, iy+1,iw,ih)], x-ix),
		    bicubicSub(p[mirror(ix-1, iy+2,iw,ih)], p[mirror(ix  , iy+2,iw,ih)], p[mirror(ix+1, iy+2,iw,ih)], p[mirror(ix+2, iy+2,iw,ih)], x-ix), y-iy);
}

float3 bicubicp(global uint32_t *p, float x, float y, int iw, int ih) {
  const int ix = (int)x, iy = (int)y;
  float3 f;
  f = bicubicSub3(bicubicSub3(unpacklab(p[mirror(ix-1, iy-1,iw,ih)]), unpacklab(p[mirror(ix  , iy-1,iw,ih)]), unpacklab(p[mirror(ix+1, iy-1,iw,ih)]), unpacklab(p[mirror(ix+2, iy-1,iw,ih)]), x-ix),
		  bicubicSub3(unpacklab(p[mirror(ix-1, iy  ,iw,ih)]), unpacklab(p[mirror(ix  , iy  ,iw,ih)]), unpacklab(p[mirror(ix+1, iy  ,iw,ih)]), unpacklab(p[mirror(ix+2, iy  ,iw,ih)]), x-ix),
		  bicubicSub3(unpacklab(p[mirror(ix-1, iy+1,iw,ih)]), unpacklab(p[mirror(ix  , iy+1,iw,ih)]), unpacklab(p[mirror(ix+1, iy+1,iw,ih)]), unpacklab(p[mirror(ix+2, iy+1,iw,ih)]), x-ix),
		  bicubicSub3(unpacklab(p[mirror(ix-1, iy+2,iw,ih)]), unpacklab(p[mirror(ix  , iy+2,iw,ih)]), unpacklab(p[mirror(ix+1, iy+2,iw,ih)]), unpacklab(p[mirror(ix+2, iy+2,iw,ih)]), x-ix), y-iy);
  return f;
}

uint32_t srgb2plab(uchar3 u) {
  const float xn = 0.950456f, zn = 1.088754f;

  int ir = s2l[u.s2];
  int ig = s2l[u.s1];
  int ib = s2l[u.s0];

  int cx = (((ir * (int)(0.412453f * 16384 + 0.5f) + ig * (int)(0.357580f * 16384 + 0.5f) + ib * (int)(0.180423f * 16384 + 0.5f)+(1<<14)) >> 15) * (int)(32768 / xn + 0.5f) + (1<<10))>>11;
  int cy = (((ir * (int)(0.212671f * 16384 + 0.5f) + ig * (int)(0.715160f * 16384 + 0.5f) + ib * (int)(0.072169f * 16384 + 0.5f)        )      )                           + (1<<10))>>11;
  int cz = (((ir * (int)(0.019334f * 16384 + 0.5f) + ig * (int)(0.119193f * 16384 + 0.5f) + ib * (int)(0.950227f * 16384 + 0.5f)+(1<<14)) >> 15) * (int)(32768 / zn + 0.5f) + (1<<10))>>11;

  int cl = (((cfunc2[cy >> 8] * (256 - (cy & 255)) + cfunc2[(cy >> 8) + 1] * (cy & 255)) >> 12) + 1) >> 1;

  int fx = cfunc[cx >> 8] * (256 - (cx & 255)) + cfunc[(cx >> 8) + 1] * (cx & 255); // <<24
  int fy = cfunc[cy >> 8] * (256 - (cy & 255)) + cfunc[(cy >> 8) + 1] * (cy & 255);
  int fz = cfunc[cz >> 8] * (256 - (cz & 255)) + cfunc[(cz >> 8) + 1] * (cz & 255);

  int fxy = (fx - fy + (1 << 7)) >> 8;
  int fyz = (fy - fz + (1 << 7)) >> 8;

  int ca = ((fxy * 8031 + (134744072 + (1 << 17))) >> 18);
  int cb = ((fyz * 3213 + (134744072 + (1 << 17))) >> 18);

  uint32_t ret = 0;
  ret = clamp(convert_uint_rtn(cb), 0U, 1023U);
  ret = (ret << 10) | clamp(convert_uint_rtn(ca), 0U, 1023U);
  ret = (ret << 12) | clamp(convert_uint_rtn(cl), 0U, 4095U);
  return ret;
}

float3 srgb2lab(uchar3 u) { return unpacklab(srgb2plab(u)); }

float icfunc(float ft) {
  if (ft > 0.20689270648f) {
    return ft * ft * ft;
  } else {
    return (ft - 16.0f/116) * (1.0f / 7.787f);
  }
}

uchar3 lab2srgb(float3 lab) {
  const float xn = 0.950456f, zn = 1.088754f;

  lab *= 256;

  float cy;
  if (lab.x > 0.20689270648f) {
    cy = (lab.x + 16) * (1.0f / 116.0f);
    cy = cy * cy * cy;
  } else {
    cy = lab.x * (1.0f / 903.3f);
  }

  float fy = (cfunc[clamp(convert_int_rtn(cy * 1024), 0, 1023)] + 9039) * (1.0f / 65536.0f);

  float fz = fy - (lab.z - 128) * (1.0f / 200.0f);
  float fx = fy + (lab.y - 128) * (1.0f / 500.0f);

  float cx = icfunc(fx) * xn;
  float cz = icfunc(fz) * zn;

  float r = cx *  3.240479f + cy * -1.537150f + cz * -0.498535f;
  float g = cx * -0.969256f + cy *  1.875991f + cz *  0.041556f;
  float b = cx *  0.055648f + cy * -0.204043f + cz *  1.057311f;

  uchar3 srgb;

  srgb.s2 = l2s[clamp(convert_int_rtn(r * 1024), 0, 1023)];
  srgb.s1 = l2s[clamp(convert_int_rtn(g * 1024), 0, 1023)];
  srgb.s0 = l2s[clamp(convert_int_rtn(b * 1024), 0, 1023)];

  return srgb;
}

uchar3 plab2srgb(uint32_t plab) { return lab2srgb(unpacklab(plab)); }

ulong xrandom(ulong s) {
  int n;
  ulong t = s;
  n = (s >> 24) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0xf3dd0fb7820fde37UL;
  n = (s >>  6) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0xe6c6ac2c59e52811UL;
  n = (s >> 18) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0x2fc7871fff7c5b45UL;
  n = (s >> 48) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0x47c7e1f70aa4f7c5UL;
  n = (s >>  0) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0x094f02b7fb9ba895UL;
  n = (s >> 12) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0x89afda817e744570UL;
  n = (s >> 36) & 63; t = ((t << n) | (t >> (64 - n))); t ^= 0xc7277d052c7bf14bUL;
  return t;
}

//

__kernel void clear(global int *out, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = 0;
}

__kernel void copy(global int *out, global int *in, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = in[x];
}

__kernel void cast_i_f(global int *out, global float *in0, float scale, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = (int)(in0[x] * scale);
}

__kernel void cast_c_i(global int8_t *out, global int *in, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = (int8_t)in[x];
}

__kernel void threshold_i_i(global int *out, global int *in, int vlow, int threshold, int vhigh, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = in[x] > threshold ? vhigh : vlow;
}

__kernel void threshold_f_f(global float *out, global float *in, float vlow, float threshold, float vhigh, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = in[x] > threshold ? vhigh : vlow;
}

__kernel void threshold_f(global float *buf, float low, float thre, float high, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  float d = buf[x];
  d = d < thre ? low : high;
  buf[x] = d;
}

__kernel void rand(global int *out, ulong seed, int size) {
  const int x = get_global_id(0);
  if (x >= size) return;

  out[x] = (int)xrandom((x    ^ 0xb21c2cb635b48285UL) * 0x9b923b9cec745401UL +
			(seed ^ 0x7bb93d75a79d2f15UL) * 0x22cab58ada573a29UL);
}

__kernel void bgr2plab(global uint32_t *out, global unsigned char *inrgb, int iw, int ih, int ws) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x, p1 = y * ws + x * 3;

  out[p0] = srgb2plab((uchar3) (inrgb[p1 + 0], inrgb[p1 + 1], inrgb[p1 + 2]));
}

__kernel void plab2bgr(global unsigned char *outrgb, global uint32_t *in, int iw, int ih, int ws) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x, p1 = y * ws + x * 3;

  uchar3 u = plab2srgb(in[p0]);
  outrgb[p1 + 0] = u.s0;
  outrgb[p1 + 1] = u.s1;
  outrgb[p1 + 2] = u.s2;
}

__kernel void convert_bgr_luminancef(global unsigned char *outrgb, global float *in, float f, int iw, int ih, int ws) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x, p1 = y * ws + x * 3;

  outrgb[p1 + 0] = outrgb[p1 + 1] = outrgb[p1 + 2] = l2s[clamp(convert_int_rtn(in[p0] * f * 1024), 0, 1023)];
}

__kernel void convert_bgr_lumaf(global unsigned char *outrgb, global float *in, float f, int iw, int ih, int ws) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x, p1 = y * ws + x * 3;

  outrgb[p1 + 0] = outrgb[p1 + 1] = outrgb[p1 + 2] = clamp(convert_int_rtn(in[p0] * f * 255), 0, 255);
}

__kernel void convert_bgr_labeli(global unsigned char *out, global int *in, int bgc, int iw, int ih, int ws) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x, p1 = y * ws + x * 3;

  int c = in[p0];

  if (c == bgc) {
    out[p1+0] = 0;
    out[p1+1] = 0;
    out[p1+2] = 0;
    return;
  }

#if 1
  int g = c * 1103515245 + 12345, d;

  d = ((g & (7 << 0)) << 5) | 31;
  out[p1+2] = (((int)d) & 255);

  d = ((g & (7 << 3)) << 2) | 31;
  out[p1+1] = (((int)d) & 255);

  d = ((g & (7 << 6)) >> 1) | 31;
  out[p1+0] = (((int)d) & 255);
#else
  out[p1+0] = ((c >>  0) & 0xff);
  out[p1+1] = ((c >>  8) & 0xff);
  out[p1+2] = ((c >> 16) & 0xff);
#endif
}

//

__kernel void pack_plab(global uint32_t *out, global float *in0, global float *in1, global float *in2, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  out[p0] = packlab((float3) (in0[p0], in1[p0], in2[p0]));
}

__kernel void unpack_plab(global float *out0, global float *out1, global float *out2, global uint32_t *in, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  float3 o = unpacklab(in[p0]);
  out0[p0] = o.x;
  out1[p0] = o.y;
  out2[p0] = o.z;
}

//

__constant float v5c[] = {
  - 4.667,-4.083, 0.000, 4.083, 4.667,
  -10.024,-0.963, 0.000, 0.963,10.024,
  -14.120, 3.622, 0.000,-3.622,14.120,
  -10.024,-0.963, 0.000, 0.963,10.024,
  - 4.667,-4.083, 0.000, 4.083, 4.667,
};

__kernel void edgevec_plab(global float2 *dst, global uint32_t *in, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  float3 vx3 = 0, vy3 = 0;

  for(int yy=-2;yy<=2;yy++) {
    for(int xx=-2;xx<=2;xx++) {
      float3 s = unpacklab(in[mirror(x+xx, y+yy,iw,ih)]);
      vx3 += v5c[(xx+2)+(yy+2)*5] * s;
      vy3 += v5c[(yy+2)+(xx+2)*5] * s;
    }
  }

  float3 ivlen3 = vx3*vx3 + vy3*vy3;

  float ivlen, vx, vy;

  if (ivlen3.s0 >= ivlen3.s1 && ivlen3.s0 >= ivlen3.s2) {
    ivlen = ivlen3.s0; vx = vx3.s0; vy = vy3.s0;
  } else if (ivlen3.s1 >= ivlen3.s2) {
    ivlen = ivlen3.s1; vx = vx3.s1; vy = vy3.s1;
  } else {
    ivlen = ivlen3.s2; vx = vx3.s2; vy = vy3.s2;
  }

  if (ivlen3.s0 >= 1e-6 && (vx3.s0 * vx + vy3.s0 * vy < 0)) {
    vx = -vx; vy = -vy;
  }

  if (ivlen > 1e-10) {
    ivlen = rsqrt(ivlen);
    vx *= ivlen; vy *= ivlen;
  } else {
    vx = vy = 0.70710678118f;
  }

  dst[p0] = (float2) (vx, vy);
}

__kernel void edgevec_f(global float2 *dst, global float *in, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  float vx = 0, vy = 0;

  for(int yy=-2;yy<=2;yy++) {
    for(int xx=-2;xx<=2;xx++) {
      float s = in[mirror(x+xx, y+yy,iw,ih)];
      vx += v5c[(xx+2)+(yy+2)*5] * s;
      vy += v5c[(yy+2)+(xx+2)*5] * s;
    }
  }

  float ivlen = vx * vx + vy * vy;

  if (ivlen > 1e-10) {
    ivlen = rsqrt(ivlen);
    vx *= ivlen; vy *= ivlen;
  } else {
    vx = vy = 0.70710678118f;
  }

  dst[p0] = (float2) (vx, vy);
}

__kernel void edge_plab(global float *out, global uint32_t *in, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;

  float3 sum = 0, t;

  t = unpacklab(in[mirror(x, y-1,iw,ih)]) + unpacklab(in[mirror(x-1, y,iw,ih)]) - unpacklab(in[mirror(x, y+1,iw,ih)]) - unpacklab(in[mirror(x+1, y,iw,ih)]);
  sum += (unpacklab(in[mirror(x-1, y-1,iw,ih)]) - unpacklab(in[mirror(x+1, y+1,iw,ih)])) * t;

  t = unpacklab(in[mirror(x, y-1,iw,ih)]) - unpacklab(in[mirror(x-1, y,iw,ih)]) + unpacklab(in[mirror(x+1, y,iw,ih)]) - unpacklab(in[mirror(x, y+1,iw,ih)]);
  sum += (unpacklab(in[mirror(x+1, y-1,iw,ih)]) - unpacklab(in[mirror(x-1, y+1,iw,ih)])) * t;

  sum = max((float3) (0, 0, 0), sum);

  out[y*iw + x] = (sum.s0 + sum.s1 + sum.s2) > 0 ? sqrt(sum.s0 + sum.s1 + sum.s2) : 0.0f;
}

__kernel void edge_f_f(global float *out, global float *in, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;

  float sum = 0, t;

  t = in[mirror(x, y-1,iw,ih)] + in[mirror(x-1, y,iw,ih)] - in[mirror(x, y+1,iw,ih)] - in[mirror(x+1, y,iw,ih)];
  sum += (in[mirror(x-1, y-1,iw,ih)] - in[mirror(x+1, y+1,iw,ih)]) * t;

  t = in[mirror(x, y-1,iw,ih)] - in[mirror(x-1, y,iw,ih)] + in[mirror(x+1, y,iw,ih)] - in[mirror(x, y+1,iw,ih)];
  sum += (in[mirror(x+1, y-1,iw,ih)] - in[mirror(x-1, y+1,iw,ih)]) * t;

  out[y*iw + x] = sqrt(max(0.0f, sum));
}

//

__kernel void thinthres_f_f_f2(global float *out, global float *in, global float2 *vxy, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  float vx = vxy[p0].x;
  float vy = vxy[p0].y;

  float am2 = bicubic(in, x-2*vx, y-2*vy, iw, ih);
  float am1 = bicubic(in, x-1*vx, y-1*vy, iw, ih);
  float a0  = in[p0];
  float ap1 = bicubic(in, x+1*vx, y+1*vy, iw, ih);
  float ap2 = bicubic(in, x+2*vx, y+2*vy, iw, ih);

  out[p0] = (am1 <= a0 && a0 >= ap1) ? (am2 + am1 + a0 + ap1 + ap2) : 0.0f;
}

__kernel void thincubic_float_float(global float *out, global float *in, global float2 *vxy, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;

  const int p0 = y * iw + x;

  float vx = vxy[p0].x;
  float vy = vxy[p0].y;

  float am2 = bicubic(in, x-2*vx, y-2*vy, iw, ih);
  float am1 = bicubic(in, x-1*vx, y-1*vy, iw, ih);
  float a0  = in[p0];
  float ap1 = bicubic(in, x+1*vx, y+1*vy, iw, ih);
  float ap2 = bicubic(in, x+2*vx, y+2*vy, iw, ih);

  const float C = 0.99f;

  out[p0] = (am2 * C <= a0 && am1 * C <= a0 && a0 >= ap1 * C && a0 >= ap2 * C) ? (am2+am1+a0+ap1+ap2) : 0;
}

//

__kernel void labelxPreprocess_int_int(global int *label, global int *pix, global int *flags, int maxpass, int bgc, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  const int p0 = y * iw + x;

  if (y == 0 && x < maxpass+1) {
    flags[x] = x == 0 ? 1 : 0;
  }

  if (x >= iw || y >= ih) return;

  if (pix[p0] == bgc) { label[p0] = -1; return; }
  if (y > 0 && pix[p0] == pix[p0-iw]) { label[p0] = p0-iw; return; }
  if (x > 0 && pix[p0] == pix[p0- 1]) { label[p0] = p0- 1; return; }
  label[p0] = p0;
}

__kernel void label8xMain_int_int(global int *label, global int *pix, global int *flags, int pass, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  if (flags[pass-1] == 0) return;

  int g = label[p0], og = g;

  if (g == -1) return;

  for(int yy=-1;yy<=1;yy++) {
    for(int xx=-1;xx<=1;xx++) {
      if (0 <=  x + xx &&  x + xx < iw && 0 <=  y + yy &&  y + yy < ih) {
	const int p1 = (y + yy) * iw + x + xx, s = label[p1];
	if (s != -1 && s < g && pix[p0] == pix[p1]) g = s;
      }
    }
  }

  for(int j=0;j<6;j++) g = label[g];

  if (g != og) {
    atomic_min(&label[og], g);
    atomic_min(&label[p0], g);
    flags[pass] = 1;
  }
}

//

__kernel void iirblur_f_f_pass0a(global float *tmp0, global float *ibuf, int r, int iw, int ih) {
  const int N = 8;
  const int y = get_global_id(0);
  if (y >= ih) return;

  constant float *coef = iircoef[r];

  float8 iv = 0.0f, tv = 0.0f;
  for(int x=-(r+1+N);x<iw;x++) {
    iv.s0 = ibuf[mirror1(x, iw) + y*iw];
    float d = iv.s0 * coef[0];
    d += coef[1] * iv.s1 + coef[2] * iv.s2 + coef[3] * iv.s3 + coef[4] * iv.s4 + coef[5] * iv.s5 + coef[6] * iv.s6 + coef[7] * iv.s7;
    d += coef[8] * tv.s0 + coef[9] * tv.s1 + coef[10] * tv.s2 + coef[11] * tv.s3 + coef[12] * tv.s4 + coef[13] * tv.s5 + coef[14] * tv.s6;
    tmp0[repeat1(x, iw) + y*iw] = d;
    iv = iv.s00123456; tv = tv.s00123456;
    tv.s0 = d;
  }
}

__kernel void iirblur_f_f_pass0b(global float *tmp1, global float *ibuf, int r, int iw, int ih) {
  const int N = 8;
  const int y = get_global_id(0);
  if (y >= ih) return;

  constant float *coef = iircoef[r];

  float8 iv = 0.0f, tv = 0.0f;
  for(int x=iw+(r+1+N);x>=0;x--) {
    iv.s0 = ibuf[mirror1(x, iw) + y*iw];
    float d = iv.s0 * coef[0];
    d += coef[1] * iv.s1 + coef[2] * iv.s2 + coef[3] * iv.s3 + coef[4] * iv.s4 + coef[5] * iv.s5 + coef[6] * iv.s6 + coef[7] * iv.s7;
    d += coef[8] * tv.s0 + coef[9] * tv.s1 + coef[10] * tv.s2 + coef[11] * tv.s3 + coef[12] * tv.s4 + coef[13] * tv.s5 + coef[14] * tv.s6;
    tmp1[repeat1(x, iw) + y*iw] = d;
    iv = iv.s00123456; tv = tv.s00123456;
    tv.s0 = d;
  }
}

__kernel void iirblur_f_f_pass1(global float *obuf, global float *tmp0, global float *tmp1, global float *ibuf, int r, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  constant float *coef = iircoef[r];

  obuf[p0] = tmp1[p0] + tmp0[p0] - ibuf[p0] * coef[0];
  tmp0[p0] = tmp1[p0] = 0;
}

__kernel void iirblur_f_f_pass2a(global float *obuf, global float *tmp0, int r, int iw, int ih) {
  const int N = 8;
  const int x = get_global_id(0);
  if (x >= iw) return;

  constant float *coef = iircoef[r];

  float8 iv = 0.0f, tv = 0.0f;
  for(int y=-(r+1+N);y<ih;y++) {
    iv.s0 = obuf[x + mirror1(y, ih)*iw];
    float d = iv.s0 * coef[0];
    d += coef[1] * iv.s1 + coef[2] * iv.s2 + coef[3] * iv.s3 + coef[4] * iv.s4 + coef[5] * iv.s5 + coef[6] * iv.s6 + coef[7] * iv.s7;
    d += coef[8] * tv.s0 + coef[9] * tv.s1 + coef[10] * tv.s2 + coef[11] * tv.s3 + coef[12] * tv.s4 + coef[13] * tv.s5 + coef[14] * tv.s6;
    tmp0[x + repeat1(y, ih)*iw] = d;
    iv = iv.s00123456; tv = tv.s00123456;
    tv.s0 = d;
  }
}

__kernel void iirblur_f_f_pass2b(global float *obuf, global float *tmp1, int r, int iw, int ih) {
  const int N = 8;
  const int x = get_global_id(0);
  if (x >= iw) return;

  constant float *coef = iircoef[r];

  float8 iv = 0.0f, tv = 0.0f;
  for(int y=ih+(r+1+N);y>=0;y--) {
    iv.s0 = obuf[x + mirror1(y, ih)*iw];
    float d = iv.s0 * coef[0];
    d += coef[1] * iv.s1 + coef[2] * iv.s2 + coef[3] * iv.s3 + coef[4] * iv.s4 + coef[5] * iv.s5 + coef[6] * iv.s6 + coef[7] * iv.s7;
    d += coef[8] * tv.s0 + coef[9] * tv.s1 + coef[10] * tv.s2 + coef[11] * tv.s3 + coef[12] * tv.s4 + coef[13] * tv.s5 + coef[14] * tv.s6;
    tmp1[x + repeat1(y, ih)*iw] = d;
    iv = iv.s00123456; tv = tv.s00123456;
    tv.s0 = d;
  }
}

__kernel void iirblur_f_f_pass3(global float *obuf, global float *tmp0, global float *tmp1, int r, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x >= iw || y >= ih) return;
  const int p0 = y * iw + x;

  constant float *coef = iircoef[r];

  obuf[p0] = tmp1[p0] + tmp0[p0] - obuf[p0] * coef[0];
}

//

__kernel void calcStrength(global int *out, global float *edge, global int *label, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x <= 0 || y <= 0 || x >= (iw-1) || y >= (ih-1)) return;
  const int p0 = y * iw + x;

  if (label[p0] <= 0) return;

  atomic_add(&out[label[p0]], (int)(edge[p0] * edge[p0] * 10000.0f));
}

__kernel void filterStrength(global int *labelinout, global int *str, int thre, int iw, int ih) {
  const int x = get_global_id(0), y = get_global_id(1);
  if (x <= 0 || y <= 0 || x >= (iw-1) || y >= (ih-1)) return;
  const int p0 = y * iw + x;

  if (labelinout[p0] <= 0 || str[labelinout[p0]] < thre) labelinout[p0] = -1;
}

//

__constant unsigned short s2l[] = {
  0,9,19,29,39,49,59,69,
  79,89,99,109,120,131,143,156,
  169,183,198,213,229,245,262,280,
  299,318,338,359,380,402,425,448,
  473,498,524,550,578,606,635,664,
  695,726,758,791,825,859,895,931,
  968,1006,1045,1084,1125,1166,1208,1251,
  1295,1340,1386,1433,1480,1529,1578,1628,
  1679,1732,1785,1839,1894,1950,2006,2064,
  2123,2183,2243,2305,2368,2431,2496,2562,
  2628,2696,2764,2834,2905,2976,3049,3123,
  3197,3273,3350,3428,3506,3586,3667,3749,
  3832,3917,4002,4088,4175,4264,4353,4444,
  4536,4628,4722,4817,4913,5011,5109,5208,
  5309,5411,5513,5617,5722,5829,5936,6044,
  6154,6265,6377,6490,6604,6720,6836,6954,
  7073,7193,7314,7437,7560,7685,7811,7939,
  8067,8197,8328,8460,8593,8727,8863,9000,
  9138,9278,9418,9560,9703,9848,9993,10140,
  10288,10438,10588,10740,10893,11048,11203,11360,
  11519,11678,11839,12001,12164,12329,12495,12662,
  12831,13000,13172,13344,13518,13693,13869,14047,
  14226,14406,14588,14771,14955,15141,15328,15516,
  15706,15897,16089,16283,16478,16675,16872,17071,
  17272,17474,17677,17882,18088,18295,18504,18714,
  18926,19139,19353,19569,19786,20004,20224,20445,
  20668,20892,21118,21345,21573,21803,22034,22267,
  22501,22736,22973,23212,23451,23693,23935,24179,
  24425,24672,24920,25170,25422,25674,25929,26184,
  26442,26700,26961,27222,27485,27750,28016,28283,
  28552,28823,29095,29368,29643,29920,30198,30477,
  30758,31041,31325,31610,31897,32186,32476,32768,
  33061,33355,33651,33949,34248
};

__constant uchar l2s[] = {
  0,3,6,10,13,15,18,20,22,24,25,27,28,30,31,32,
  34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,
  50,50,51,52,53,54,54,55,56,57,57,58,59,59,60,61,
  61,62,63,63,64,65,65,66,66,67,68,68,69,69,70,70,
  71,72,72,73,73,74,74,75,75,76,76,77,77,78,78,79,
  79,80,80,81,81,82,82,83,83,84,84,84,85,85,86,86,
  87,87,88,88,88,89,89,90,90,91,91,91,92,92,93,93,
  93,94,94,95,95,95,96,96,97,97,97,98,98,98,99,99,
  100,100,100,101,101,101,102,102,102,103,103,103,104,104,105,105,
  105,106,106,106,107,107,107,108,108,108,109,109,109,110,110,110,
  111,111,111,112,112,112,113,113,113,113,114,114,114,115,115,115,
  116,116,116,117,117,117,117,118,118,118,119,119,119,120,120,120,
  120,121,121,121,122,122,122,122,123,123,123,124,124,124,124,125,
  125,125,126,126,126,126,127,127,127,127,128,128,128,129,129,129,
  129,130,130,130,130,131,131,131,131,132,132,132,133,133,133,133,
  134,134,134,134,135,135,135,135,136,136,136,136,137,137,137,137,
  138,138,138,138,139,139,139,139,140,140,140,140,140,141,141,141,
  141,142,142,142,142,143,143,143,143,144,144,144,144,144,145,145,
  145,145,146,146,146,146,147,147,147,147,147,148,148,148,148,149,
  149,149,149,149,150,150,150,150,151,151,151,151,151,152,152,152,
  152,153,153,153,153,153,154,154,154,154,154,155,155,155,155,156,
  156,156,156,156,157,157,157,157,157,158,158,158,158,158,159,159,
  159,159,159,160,160,160,160,161,161,161,161,161,162,162,162,162,
  162,163,163,163,163,163,164,164,164,164,164,164,165,165,165,165,
  165,166,166,166,166,166,167,167,167,167,167,168,168,168,168,168,
  169,169,169,169,169,169,170,170,170,170,170,171,171,171,171,171,
  172,172,172,172,172,172,173,173,173,173,173,174,174,174,174,174,
  175,175,175,175,175,175,176,176,176,176,176,176,177,177,177,177,
  177,178,178,178,178,178,178,179,179,179,179,179,179,180,180,180,
  180,180,181,181,181,181,181,181,182,182,182,182,182,182,183,183,
  183,183,183,183,184,184,184,184,184,184,185,185,185,185,185,186,
  186,186,186,186,186,187,187,187,187,187,187,188,188,188,188,188,
  188,188,189,189,189,189,189,189,190,190,190,190,190,190,191,191,
  191,191,191,191,192,192,192,192,192,192,193,193,193,193,193,193,
  194,194,194,194,194,194,194,195,195,195,195,195,195,196,196,196,
  196,196,196,197,197,197,197,197,197,197,198,198,198,198,198,198,
  199,199,199,199,199,199,199,200,200,200,200,200,200,201,201,201,
  201,201,201,201,202,202,202,202,202,202,202,203,203,203,203,203,
  203,204,204,204,204,204,204,204,205,205,205,205,205,205,205,206,
  206,206,206,206,206,206,207,207,207,207,207,207,207,208,208,208,
  208,208,208,208,209,209,209,209,209,209,209,210,210,210,210,210,
  210,210,211,211,211,211,211,211,211,212,212,212,212,212,212,212,
  213,213,213,213,213,213,213,214,214,214,214,214,214,214,215,215,
  215,215,215,215,215,216,216,216,216,216,216,216,217,217,217,217,
  217,217,217,217,218,218,218,218,218,218,218,219,219,219,219,219,
  219,219,219,220,220,220,220,220,220,220,221,221,221,221,221,221,
  221,222,222,222,222,222,222,222,222,223,223,223,223,223,223,223,
  223,224,224,224,224,224,224,224,225,225,225,225,225,225,225,225,
  226,226,226,226,226,226,226,226,227,227,227,227,227,227,227,228,
  228,228,228,228,228,228,228,229,229,229,229,229,229,229,229,230,
  230,230,230,230,230,230,230,231,231,231,231,231,231,231,231,232,
  232,232,232,232,232,232,232,233,233,233,233,233,233,233,233,234,
  234,234,234,234,234,234,234,235,235,235,235,235,235,235,235,236,
  236,236,236,236,236,236,236,237,237,237,237,237,237,237,237,238,
  238,238,238,238,238,238,238,238,239,239,239,239,239,239,239,239,
  240,240,240,240,240,240,240,240,241,241,241,241,241,241,241,241,
  241,242,242,242,242,242,242,242,242,243,243,243,243,243,243,243,
  243,243,244,244,244,244,244,244,244,244,245,245,245,245,245,245,
  245,245,245,246,246,246,246,246,246,246,246,247,247,247,247,247,
  247,247,247,247,248,248,248,248,248,248,248,248,248,249,249,249,
  249,249,249,249,249,250,250,250,250,250,250,250,250,250,251,251,
  251,251,251,251,251,251,251,252,252,252,252,252,252,252,252,252,
  253,253,253,253,253,253,253,253,253,254,254,254,254,254,254,254,
  254,254,255,255,255,255,255,255,255,255,255,255,255,255,255,255,
};

__constant unsigned short cfunc[] = {
  0,499,997,1496,1994,2492,2991,3489,3987,4486,4969,5421,5847,6249,6632,6996,
  7345,7679,8001,8311,8610,8900,9180,9452,9716,9973,10223,10467,10705,10937,11164,11386,
  11604,11816,12025,12229,12430,12627,12821,13011,13198,13381,13562,13740,13915,14088,14258,14426,
  14591,14754,14915,15073,15230,15384,15537,15688,15837,15984,16129,16273,16415,16556,16695,16833,
  16969,17104,17237,17369,17500,17629,17758,17885,18010,18135,18259,18381,18502,18623,18742,18860,
  18977,19093,19209,19323,19437,19549,19661,19772,19882,19991,20099,20207,20313,20419,20525,20629,
  20733,20836,20938,21040,21141,21241,21340,21439,21538,21635,21732,21829,21925,22020,22115,22209,
  22302,22395,22488,22580,22671,22762,22852,22942,23032,23120,23209,23297,23384,23471,23557,23643,
  23729,23814,23899,23983,24067,24150,24233,24316,24398,24480,24561,24642,24723,24803,24883,24962,
  25041,25120,25198,25276,25354,25431,25508,25585,25661,25737,25812,25888,25963,26037,26112,26186,
  26259,26333,26406,26479,26551,26623,26695,26767,26838,26909,26980,27050,27121,27190,27260,27330,
  27399,27468,27536,27605,27673,27741,27808,27875,27943,28009,28076,28143,28209,28275,28340,28406,
  28471,28536,28601,28665,28730,28794,28858,28921,28985,29048,29111,29174,29237,29299,29361,29423,
  29485,29547,29608,29670,29731,29792,29852,29913,29973,30033,30093,30153,30212,30272,30331,30390,
  30449,30507,30566,30624,30682,30740,30798,30856,30913,30971,31028,31085,31142,31198,31255,31311,
  31367,31423,31479,31535,31591,31646,31701,31757,31811,31866,31921,31976,32030,32084,32138,32192,
  32246,32300,32353,32407,32460,32513,32566,32619,32672,32724,32777,32829,32881,32933,32985,33037,
  33089,33140,33192,33243,33294,33345,33396,33447,33498,33548,33599,33649,33700,33750,33800,33849,
  33899,33949,33998,34048,34097,34146,34195,34244,34293,34342,34391,34439,34487,34536,34584,34632,
  34680,34728,34776,34823,34871,34918,34966,35013,35060,35107,35154,35201,35248,35295,35341,35388,
  35434,35480,35526,35573,35619,35664,35710,35756,35802,35847,35893,35938,35983,36028,36073,36118,
  36163,36208,36253,36297,36342,36386,36431,36475,36519,36563,36607,36651,36695,36739,36782,36826,
  36870,36913,36956,37000,37043,37086,37129,37172,37215,37258,37300,37343,37385,37428,37470,37513,
  37555,37597,37639,37681,37723,37765,37807,37848,37890,37932,37973,38015,38056,38097,38138,38180,
  38221,38262,38303,38343,38384,38425,38466,38506,38547,38587,38627,38668,38708,38748,38788,38828,
  38868,38908,38948,38988,39027,39067,39106,39146,39185,39225,39264,39303,39342,39382,39421,39460,
  39499,39537,39576,39615,39654,39692,39731,39769,39808,39846,39884,39923,39961,39999,40037,40075,
  40113,40151,40189,40227,40264,40302,40340,40377,40415,40452,40489,40527,40564,40601,40638,40675,
  40712,40749,40786,40823,40860,40897,40934,40970,41007,41043,41080,41116,41153,41189,41225,41262,
  41298,41334,41370,41406,41442,41478,41514,41550,41585,41621,41657,41693,41728,41764,41799,41835,
  41870,41905,41941,41976,42011,42046,42081,42116,42151,42186,42221,42256,42291,42325,42360,42395,
  42429,42464,42498,42533,42567,42602,42636,42670,42705,42739,42773,42807,42841,42875,42909,42943,
  42977,43011,43045,43078,43112,43146,43179,43213,43246,43280,43313,43347,43380,43414,43447,43480,
  43513,43546,43580,43613,43646,43679,43712,43744,43777,43810,43843,43876,43908,43941,43974,44006,
  44039,44071,44104,44136,44169,44201,44233,44265,44298,44330,44362,44394,44426,44458,44490,44522,
  44554,44586,44618,44650,44681,44713,44745,44777,44808,44840,44871,44903,44934,44966,44997,45028,
  45060,45091,45122,45154,45185,45216,45247,45278,45309,45340,45371,45402,45433,45464,45495,45525,
  45556,45587,45618,45648,45679,45709,45740,45770,45801,45831,45862,45892,45923,45953,45983,46013,
  46044,46074,46104,46134,46164,46194,46224,46254,46284,46314,46344,46374,46404,46433,46463,46493,
  46523,46552,46582,46612,46641,46671,46700,46730,46759,46788,46818,46847,46877,46906,46935,46964,
  46993,47023,47052,47081,47110,47139,47168,47197,47226,47255,47284,47313,47342,47370,47399,47428,
  47457,47485,47514,47543,47571,47600,47628,47657,47685,47714,47742,47771,47799,47827,47856,47884,
  47912,47940,47969,47997,48025,48053,48081,48109,48137,48165,48193,48221,48249,48277,48305,48333,
  48361,48388,48416,48444,48472,48499,48527,48555,48582,48610,48637,48665,48692,48720,48747,48775,
  48802,48830,48857,48884,48912,48939,48966,48993,49020,49048,49075,49102,49129,49156,49183,49210,
  49237,49264,49291,49318,49345,49372,49399,49425,49452,49479,49506,49532,49559,49586,49612,49639,
  49666,49692,49719,49745,49772,49798,49825,49851,49878,49904,49930,49957,49983,50009,50036,50062,
  50088,50114,50140,50167,50193,50219,50245,50271,50297,50323,50349,50375,50401,50427,50453,50479,
  50504,50530,50556,50582,50608,50633,50659,50685,50710,50736,50762,50787,50813,50838,50864,50890,
  50915,50941,50966,50991,51017,51042,51068,51093,51118,51144,51169,51194,51219,51245,51270,51295,
  51320,51345,51370,51396,51421,51446,51471,51496,51521,51546,51571,51596,51620,51645,51670,51695,
  51720,51745,51770,51794,51819,51844,51868,51893,51918,51942,51967,51992,52016,52041,52065,52090,
  52114,52139,52163,52188,52212,52237,52261,52285,52310,52334,52358,52383,52407,52431,52456,52480,
  52504,52528,52552,52576,52601,52625,52649,52673,52697,52721,52745,52769,52793,52817,52841,52865,
  52889,52913,52936,52960,52984,53008,53032,53055,53079,53103,53127,53150,53174,53198,53221,53245,
  53269,53292,53316,53339,53363,53386,53410,53433,53457,53480,53504,53527,53551,53574,53597,53621,
  53644,53667,53691,53714,53737,53760,53784,53807,53830,53853,53876,53899,53923,53946,53969,53992,
  54015,54038,54061,54084,54107,54130,54153,54176,54199,54222,54244,54267,54290,54313,54336,54359,
  54381,54404,54427,54450,54472,54495,54518,54541,54563,54586,54608,54631,54654,54676,54699,54721,
  54744,54766,54789,54811,54834,54856,54879,54901,54924,54946,54968,54991,55013,55035,55058,55080,
  55102,55124,55147,55169,55191,55213,55236,55258,55280,55302,55324,55346,55368,55390,55412,55435,
  55457,55479,55501,55523,55545,55567,55588,55610,55632,55654,55676,55698,55720,55742,55763,55785,
  55807,55829,55851,55872,55894,55916,55938,55959,55981,56003,56024,56046,56068,56089,56111,56132,
  56154,56175,56197,56218,56240,56261,56283,56304,56326,56347,56369,56390,56412,56433,56454,56476,
  56497,56518,56540,56561,56582,56603,56625
};

__constant unsigned short cfunc2[] = {
  0,453,907,1360,1814,2267,2721,3174,3627,4081,4521,4932,5319,5685,6033,6365,
  6682,6986,7279,7561,7833,8096,8351,8599,8839,9073,9301,9522,9739,9950,10157,10359,
  10557,10750,10940,11126,11309,11488,11664,11837,12007,12174,12338,12500,12660,12817,12972,13124,
  13274,13423,13569,13713,13856,13996,14135,14272,14408,14542,14674,14805,14934,15062,15189,15314,
  15438,15561,15682,15802,15921,16039,16156,16271,16386,16499,16611,16723,16833,16942,17051,17159,
  17265,17371,17476,17580,17683,17785,17887,17988,18088,18187,18286,18384,18481,18577,18673,18768,
  18862,18956,19049,19142,19233,19325,19415,19505,19595,19684,19772,19860,19947,20034,20120,20205,
  20290,20375,20459,20543,20626,20709,20791,20873,20954,21035,21115,21195,21274,21354,21432,21510,
  21588,21666,21743,21819,21896,21972,22047,22122,22197,22271,22345,22419,22492,22565,22638,22710,
  22782,22854,22925,22996,23067,23137,23207,23277,23346,23415,23484,23552,23620,23688,23756,23823,
  23890,23957,24024,24090,24156,24221,24287,24352,24417,24482,24546,24610,24674,24738,24801,24864,
  24927,24990,25052,25114,25176,25238,25300,25361,25422,25483,25543,25604,25664,25724,25784,25843,
  25903,25962,26021,26079,26138,26196,26254,26312,26370,26428,26485,26542,26599,26656,26713,26769,
  26825,26881,26937,26993,27049,27104,27159,27214,27269,27324,27378,27433,27487,27541,27595,27648,
  27702,27755,27809,27862,27915,27967,28020,28072,28125,28177,28229,28281,28332,28384,28435,28487,
  28538,28589,28640,28690,28741,28791,28842,28892,28942,28992,29041,29091,29141,29190,29239,29288,
  29337,29386,29435,29483,29532,29580,29628,29676,29724,29772,29820,29868,29915,29963,30010,30057,
  30104,30151,30198,30244,30291,30337,30384,30430,30476,30522,30568,30614,30660,30705,30751,30796,
  30841,30886,30931,30976,31021,31066,31111,31155,31200,31244,31288,31332,31376,31420,31464,31508,
  31552,31595,31639,31682,31725,31769,31812,31855,31898,31940,31983,32026,32068,32111,32153,32195,
  32238,32280,32322,32364,32405,32447,32489,32531,32572,32613,32655,32696,32737,32778,32819,32860,
  32901,32942,32982,33023,33064,33104,33144,33185,33225,33265,33305,33345,33385,33425,33464,33504,
  33544,33583,33623,33662,33701,33741,33780,33819,33858,33897,33936,33974,34013,34052,34090,34129,
  34167,34206,34244,34282,34320,34358,34396,34434,34472,34510,34548,34585,34623,34661,34698,34736,
  34773,34810,34847,34885,34922,34959,34996,35033,35069,35106,35143,35180,35216,35253,35289,35326,
  35362,35398,35434,35471,35507,35543,35579,35615,35651,35686,35722,35758,35794,35829,35865,35900,
  35936,35971,36006,36041,36077,36112,36147,36182,36217,36252,36287,36321,36356,36391,36425,36460,
  36495,36529,36563,36598,36632,36666,36701,36735,36769,36803,36837,36871,36905,36939,36973,37006,
  37040,37074,37107,37141,37174,37208,37241,37274,37308,37341,37374,37407,37441,37474,37507,37540,
  37573,37605,37638,37671,37704,37736,37769,37802,37834,37867,37899,37932,37964,37996,38029,38061,
  38093,38125,38157,38189,38221,38253,38285,38317,38349,38381,38412,38444,38476,38507,38539,38571,
  38602,38633,38665,38696,38728,38759,38790,38821,38852,38884,38915,38946,38977,39008,39038,39069,
  39100,39131,39162,39192,39223,39254,39284,39315,39345,39376,39406,39437,39467,39497,39528,39558,
  39588,39618,39648,39678,39709,39739,39769,39798,39828,39858,39888,39918,39948,39977,40007,40037,
  40066,40096,40125,40155,40184,40214,40243,40273,40302,40331,40360,40390,40419,40448,40477,40506,
  40535,40564,40593,40622,40651,40680,40709,40737,40766,40795,40824,40852,40881,40910,40938,40967,
  40995,41024,41052,41080,41109,41137,41165,41194,41222,41250,41278,41306,41335,41363,41391,41419,
  41447,41475,41503,41530,41558,41586,41614,41642,41669,41697,41725,41752,41780,41808,41835,41863,
  41890,41918,41945,41973,42000,42027,42055,42082,42109,42136,42164,42191,42218,42245,42272,42299,
  42326,42353,42380,42407,42434,42461,42488,42514,42541,42568,42595,42621,42648,42675,42701,42728,
  42754,42781,42808,42834,42860,42887,42913,42940,42966,42992,43019,43045,43071,43097,43123,43150,
  43176,43202,43228,43254,43280,43306,43332,43358,43384,43410,43436,43461,43487,43513,43539,43565,
  43590,43616,43642,43667,43693,43719,43744,43770,43795,43821,43846,43872,43897,43922,43948,43973,
  43998,44024,44049,44074,44099,44125,44150,44175,44200,44225,44250,44275,44300,44325,44350,44375,
  44400,44425,44450,44475,44500,44524,44549,44574,44599,44623,44648,44673,44697,44722,44747,44771,
  44796,44820,44845,44869,44894,44918,44943,44967,44991,45016,45040,45064,45089,45113,45137,45161,
  45186,45210,45234,45258,45282,45306,45330,45354,45378,45402,45426,45450,45474,45498,45522,45546,
  45570,45594,45618,45641,45665,45689,45713,45736,45760,45784,45807,45831,45854,45878,45902,45925,
  45949,45972,45996,46019,46043,46066,46089,46113,46136,46159,46183,46206,46229,46253,46276,46299,
  46322,46345,46369,46392,46415,46438,46461,46484,46507,46530,46553,46576,46599,46622,46645,46668,
  46691,46714,46737,46759,46782,46805,46828,46851,46873,46896,46919,46941,46964,46987,47009,47032,
  47055,47077,47100,47122,47145,47167,47190,47212,47235,47257,47279,47302,47324,47347,47369,47391,
  47414,47436,47458,47480,47503,47525,47547,47569,47591,47613,47636,47658,47680,47702,47724,47746,
  47768,47790,47812,47834,47856,47878,47900,47922,47943,47965,47987,48009,48031,48053,48074,48096,
  48118,48140,48161,48183,48205,48226,48248,48270,48291,48313,48334,48356,48378,48399,48421,48442,
  48464,48485,48506,48528,48549,48571,48592,48613,48635,48656,48677,48699,48720,48741,48763,48784,
  48805,48826,48847,48869,48890,48911,48932,48953,48974,48995,49016,49037,49059,49080,49101,49122,
  49143,49163,49184,49205,49226,49247,49268,49289,49310,49331,49351,49372,49393,49414,49435,49455,
  49476,49497,49517,49538,49559,49579,49600,49621,49641,49662,49683,49703,49724,49744,49765,49785,
  49806,49826,49847,49867,49888,49908,49928,49949,49969,49990,50010,50030,50051,50071,50091,50112,
  50132,50152,50172,50193,50213,50233,50253,50273,50293,50314,50334,50354,50374,50394,50414,50434,
  50454,50474,50494,50514,50534,50554,50574,50594,50614,50634,50654,50674,50694,50714,50733,50753,
  50773,50793,50813,50833,50852,50872,50892,50912,50931,50951,50971,50990,51010,51030,51049,51069,
  51089,51108,51128,51147,51167,51187,51206,51226,51245,51265,51284,51304,51323,51342,51362,51381,
  51401,51420,51440,51459,51478,51498,51517
};

__constant float iircoef[32][15] = {
  // 1
  {
    1.1968268412f, 0.0066482081f, -0.0000369159f, 0.0000002050f, 
    -0.0000000011f, 0.0000000000f, -0.0000000000f, 0.0000000000f, 
    0.0055541344f, -0.0000308409f, 0.0000001713f, -0.0000000010f, 
    0.0000000000f, -0.0000000000f, -0.0000000000f, 
  },
  // 2
  {
    0.5984134206f, 0.1038383077f, -0.0121238572f, 0.0013108988f, 
    -0.0001376963f, 0.0000144192f, -0.0000023625f, -0.0000000267f, 
    0.1511297752f, -0.0176956557f, 0.0019154795f, -0.0002012217f, 
    0.0000206591f, -0.0000006002f, 0.0000000180f, 
  },
  // 3
  {
    0.3989422804f, 0.1414542400f, -0.0030406818f, -0.0041116157f, 
    0.0006696623f, 0.0000498707f, -0.0000449761f, -0.0000051528f, 
    0.2519574622f, -0.0098627835f, -0.0067013653f, 0.0012572396f, 
    0.0000481394f, -0.0000097781f, 0.0000006462f, 
  },
  // 4
  {
    0.2992067105f, 0.1366767807f, 0.0208926040f, -0.0069574878f, 
    -0.0020213385f, 0.0005085353f, 0.0003107938f, 0.0000374566f, 
    0.2980424247f, 0.0298515862f, -0.0164807453f, -0.0030985243f, 
    0.0011876876f, -0.0001730692f, 0.0000121781f, 
  },
  // 5
  {
    0.2393653696f, 0.1233259528f, 0.0361222057f, -0.0008004665f, 
    -0.0051378852f, -0.0028052364f, -0.0004582368f, -0.0000317911f, 
    0.3200496611f, 0.0685160638f, -0.0117714980f, -0.0092560053f, 
    0.0047643905f, -0.0010979703f, 0.0001190780f, 
  },
  // 6
  {
    0.1994711473f, 0.1097468905f, 0.0425939826f, 0.0067663095f, 
    -0.0033515935f, -0.0026165692f, -0.0016459693f, -0.0002414118f, 
    0.3323078887f, 0.0997363125f, 0.0011591449f, -0.0172636119f, 
    -0.0057667710f, 0.0064460475f, -0.0014183079f, 
  },
  // 7
  {
    0.1709750735f, 0.0984663399f, 0.0459142899f, 0.0167696504f, 
    0.0065166958f, 0.0028273657f, -0.0010413533f, -0.0001359209f, 
    0.3363447623f, 0.1171937469f, -0.0003709590f, -0.0360471291f, 
    -0.0113865917f, 0.0174437637f, -0.0043934342f, 
  },
  // 8
  {
    0.1496033042f, 0.0889772544f, 0.0466966229f, 0.0226458703f, 
    0.0126167248f, 0.0073512509f, 0.0002205437f, 0.0001354835f, 
    0.3373501960f, 0.1282570687f, 0.0055357179f, -0.0408180282f, 
    -0.0204814964f, 0.0252468306f, -0.0061513566f, 
  },
  // 9
  {
    0.1329803941f, 0.0808595083f, 0.0443823940f, 0.0225982826f, 
    0.0128388121f, 0.0083652261f, 0.0019791649f, 0.0003865980f, 
    0.3379115748f, 0.1473327968f, 0.0266207466f, -0.0335193079f, 
    -0.0314352674f, 0.0160182304f, -0.0002276320f, 
  },
  // 10
  {
    0.1196820170f, 0.0723369053f, 0.0388804292f, 0.0186848390f, 
    0.0101468368f, 0.0077015455f, 0.0036393790f, 0.0006022105f, 
    0.3515969983f, 0.1742863583f, 0.0505556059f, -0.0264340529f, 
    -0.0440305905f, -0.0010337684f, 0.0098366617f, 
  },
  // 11
  {
    0.1088004131f, 0.0573636997f, 0.0226045985f, 0.0063100547f, 
    0.0057236804f, 0.0090285503f, 0.0024463416f, 0.0005349718f, 
    0.4363356202f, 0.2334577615f, 0.0565680941f, -0.0688967536f, 
    -0.0785915265f, 0.0334050368f, 0.0038539731f, 
  },
  // 12
  {
    0.0997357385f, 0.0254438064f, -0.0144607797f, -0.0156748784f, 
    0.0076179980f, 0.0194933891f, -0.0068375272f, 0.0001311987f, 
    0.7141211197f, 0.3353369335f, -0.0432259066f, -0.2629365938f, 
    -0.1308849355f, 0.2543954490f, -0.0767902101f, 
  },
  // 13
  {
    0.0920630989f, 0.0327952390f, 0.0001272727f, -0.0065098551f, 
    0.0033483906f, 0.0105441981f, 0.0004310366f, 0.0004636846f, 
    0.6174921445f, 0.2963384295f, 0.0139684376f, -0.1491830856f, 
    -0.1043796734f, 0.0926402634f, -0.0109533693f, 
  },
  // 14
  {
    0.0854868549f, 0.0412955715f, 0.0104850922f, -0.0063088790f, 
    -0.0086460350f, 0.0004000430f, 0.0104696817f, 0.0004660982f, 
    0.4942425427f, 0.3065505734f, 0.1366446441f, -0.0214618571f, 
    -0.1367013973f, -0.1335800264f, 0.1083488761f, 
  },
  // 15
  {
    0.0797886052f, -0.0045040261f, -0.0323961412f, -0.0165407784f, 
    0.0136018434f, 0.0201893627f, -0.0101715888f, 0.0000663388f, 
    1.0366487839f, 0.3130194448f, -0.2211908975f, -0.3823469815f, 
    -0.0817610185f, 0.3760241086f, -0.1330920811f, 
  },
  // 16
  {
    0.0748018489f, 0.0261185380f, -0.0024399774f, -0.0124734514f, 
    -0.0081594062f, 0.0025033765f, 0.0085728400f, 0.0003713940f, 
    0.6334164040f, 0.3423828568f, 0.0936442351f, -0.0880049759f, 
    -0.1603219845f, -0.0905383247f, 0.1034372150f, 
  },
  // 17
  {
    0.0704014767f, -0.0148022032f, -0.0355572407f, -0.0137463850f, 
    0.0163541684f, 0.0194032960f, -0.0112323176f, 0.0000440127f, 
    1.1948031041f, 0.2683341589f, -0.3223485288f, -0.4261529631f, 
    -0.0401726093f, 0.4319095505f, -0.1640627566f, 
  },
  // 18
  {
    0.0664903264f, -0.0191764187f, -0.0367504016f, -0.0121423990f, 
    0.0181146894f, 0.0193702480f, -0.0121766674f, 0.0000315479f, 
    1.2746150096f, 0.2416435676f, -0.3789271494f, -0.4514287373f, 
    -0.0149022723f, 0.4716010196f, -0.1871837758f, 
  },
  // 19
  {
    0.0629909055f, -0.0139764662f, -0.0311630337f, -0.0121449026f, 
    0.0125759666f, 0.0153094387f, -0.0075397338f, 0.0000618978f, 
    1.2095015170f, 0.2515757382f, -0.3124437293f, -0.3923665185f, 
    -0.0417175759f, 0.3601516008f, -0.1237790604f, 
  },
  // 20
  {
    0.0598413309f, -0.0240888161f, -0.0364018545f, -0.0091189615f, 
    0.0196346378f, 0.0182935803f, -0.0127143078f, 0.0000206271f, 
    1.3913565154f, 0.1885126880f, -0.4604429309f, -0.4751463221f, 
    0.0266317987f, 0.5157431561f, -0.2158616818f, 
  },
  // 21
  {
    0.0569924804f, 0.0113490242f, -0.0130001782f, -0.0190479222f, 
    -0.0114353393f, 0.0023891792f, 0.0113210030f, 0.0001482904f, 
    0.7907048940f, 0.4054153690f, 0.0860642538f, -0.1457230772f, 
    -0.2469201043f, -0.1573042787f, 0.1947659290f, 
  },
  // 22
  {
    0.0544018588f, 0.0038413128f, -0.0172152173f, -0.0172641130f, 
    -0.0063507678f, 0.0048320006f, 0.0076071097f, 0.0001092293f, 
    0.9201526804f, 0.3683143414f, -0.0143859690f, -0.2083920960f, 
    -0.2076264668f, -0.0513000272f, 0.1366425702f, 
  },
  // 23
  {
    0.0520358718f, 0.0089797805f, -0.0139243775f, -0.0194274320f, 
    -0.0116981619f, 0.0023618595f, 0.0117836873f, 0.0001256072f, 
    0.8189491307f, 0.4221318247f, 0.0895417511f, -0.1577994690f, 
    -0.2728036834f, -0.1806756884f, 0.2234128600f, 
  },
  // 24
  {
    0.0498677919f, -0.0285695679f, -0.0336428417f, -0.0043820623f, 
    0.0209224946f, 0.0161412805f, -0.0130880847f, 0.0000093660f, 
    1.5651241194f, 0.0909293622f, -0.5872168740f, -0.5014047951f, 
    0.0995763881f, 0.5841367859f, -0.2649874889f, 
  },
  // 25
  {
    0.0478735114f, -0.0041376428f, -0.0198738559f, -0.0138254816f, 
    -0.0008582127f, 0.0065113819f, 0.0034108657f, 0.0000862718f, 
    1.0791767807f, 0.3153115446f, -0.1355874884f, -0.2738326725f, 
    -0.1546621584f, 0.0647741526f, 0.0682652094f, 
  },
  // 26
  {
    0.0460318422f, -0.0294730777f, -0.0320325317f, -0.0025726498f, 
    0.0211383349f, 0.0151784312f, -0.0131393357f, 0.0000063673f, 
    1.6336425739f, 0.0468000939f, -0.6394697088f, -0.5092315092f, 
    0.1328040841f, 0.6132983924f, -0.2876784450f, 
  },
  // 27
  {
    0.0443269440f, -0.0285176779f, -0.0307368629f, -0.0024221005f, 
    0.0201074494f, 0.0143707157f, -0.0123189581f, 0.0000073018f, 
    1.6371952610f, 0.0419056540f, -0.6383107865f, -0.5028818447f, 
    0.1324628924f, 0.6004006350f, -0.2799989041f, 
  },
  // 28
  {
    0.0427433279f, 0.0032390349f, -0.0159814229f, -0.0187791797f, 
    -0.0100542794f, 0.0033000240f, 0.0113682253f, 0.0000592097f, 
    0.9185005748f, 0.4379588120f, 0.0559255865f, -0.2083625246f, 
    -0.3121557527f, -0.1869821969f, 0.2648587870f, 
  },
  // 29
  {
    0.0412687184f, -0.0020754079f, -0.0182012088f, -0.0162521312f, 
    -0.0055456600f, 0.0050013013f, 0.0080916593f, 0.0000419546f, 
    1.0450223591f, 0.3804711425f, -0.0545434003f, -0.2617866513f, 
    -0.2545069277f, -0.0733756035f, 0.1952470873f, 
  },
  // 30
  {
    0.0398937661f, -0.0022548686f, -0.0177970661f, -0.0157912041f, 
    -0.0053646477f, 0.0048462754f, 0.0078356262f, 0.0000388455f, 
    1.0515598741f, 0.3800115152f, -0.0569951170f, -0.2634621868f, 
    -0.2549635212f, -0.0735128318f, 0.1955608148f, 
  },
  // 31
  {
    0.0386073195f, -0.0026573620f, -0.0177563852f, -0.0156267771f, 
    -0.0052464524f, 0.0049054277f, 0.0079060938f, 0.0000345602f, 
    1.0641533107f, 0.3821771594f, -0.0613065680f, -0.2705980722f, 
    -0.2617739706f, -0.0761864939f, 0.2041190380f, 
  },
  // 32
  {
    0.0374010952f, 0.0011188355f, -0.0159050509f, -0.0176196657f, 
    -0.0089191574f, 0.0035997631f, 0.0107972030f, 0.0000372860f, 
    0.9657127782f, 0.4463648375f, 0.0390300662f, -0.2351516544f, 
    -0.3337316692f, -0.1905108661f, 0.2881789365f, 
  },
};