void filterExpand(Img *dst, Img *src, int l, int r, int t, int b) {

  if (!src->data) return;

  if (l < 0) l = 0;

  if (r < 0) r = 0;

  if (t < 0) t = 0;

  if (b < 0) b = 0;

  int w = src->w + l + r;

  int h = src->h + t + b;

  if (dst->w != w || dst->h != h) imgInit(dst, w, h);

  double *ps = src->data, *pd = dst->data + w*t*4;

  for(int y = 0; y < src->h; ++y) {

    memcpy(pd, ps, 4*sizeof(double)*src->w);

    for(double *p = ps, *e = p + 4*l; p < e; p += 4)

      memcpy(p, pd, 4*sizeof(double));

    ps += (src->w-1)*4;

    pd += w*4;

    for(double *e = pd, *p = e - r*4; p < e; p += 4)

      memcpy(p, ps, 4*sizeof(double));

    ps += 4;

  }

  for(int y = 0; y < t; ++y)

    memcpy(imgPixel(dst, 0, y), imgPixel(dst, 0, t), 4*sizeof(double)*w);

  for(int y = h-b; y < h; ++y)

    memcpy(imgPixel(dst, 0, y), imgPixel(dst, 0, h-b-1), 4*sizeof(double)*w);

}

void filterArea(Img *dst, Img *src) {

  int w = dst->w > src->w ? dst->w : src->w;

  int h = dst->h > src->h ? dst->h : src->h;

  if (dst->w < w || dst->h < h)

    imgInit(dst, w, h);

  int gap = (w - src->w)*4;

  double *pd = dst->data, *ps = src->data;

  memcpy(pd, ps, 4*sizeof(*ps)); ps += 4; pd += 4;

  for(double *e = ps + 4*(src->w - 1); ps < e; ++ps, ++pd)

    *pd = *ps + *(pd - 4);

  pd += gap;

  for(int y = 1; y < src->h; ++y) {

    for(int i = 0; i < 4; ++i, ++ps, ++pd)

      *pd = *ps + *(pd - w*4);

    for(double *e = ps + 4*(src->w - 1); ps < e; ++ps, ++pd)

      *pd = *ps + *(pd - 4) + *(pd - w*4) - *(pd - (w+1)*4);

    pd += gap;

  }

}

void pixelAreaAvg(double *pixel, Img *areaImg, int x0, int y0, int x1, int y1) {

  if (x0 < 0) x0 = 0;

  if (y0 < 0) y0 = 0;

  if (x1 >= areaImg->w) x1 = areaImg->w - 1;

  if (y1 >= areaImg->h) y1 = areaImg->h - 1;

  if (x1 < x0 || y1 < y0) { memset(pixel, 0, sizeof(*pixel)*4); return; }

  double k = 1.0/((x1 - x0 + 1)*(y1 - y0 + 1));

  if (x0) {

    if (y0) {

      --x0; --y0;

      double *p00 = imgPixel(areaImg, x0, y0);

      double *p01 = imgPixel(areaImg, x0, y1);

      double *p10 = imgPixel(areaImg, x1, y0);

      double *p11 = imgPixel(areaImg, x1, y1);

      for(int i = 0; i < 4; ++i)

        pixel[i] = (p00[i] + p11[i] - p01[i] - p10[i])*k;

      return;

    }

    double *p0 = imgPixel(areaImg, x0-1, y1);

    double *p1 = imgPixel(areaImg, x1, y1);

    for(int i = 0; i < 4; ++i)

      pixel[i] = (p1[i] - p0[i])*k;

  } else

  if (y0) {

    double *p0 = imgPixel(areaImg, x1, y0-1);

    double *p1 = imgPixel(areaImg, x1, y1);

    for(int i = 0; i < 4; ++i)

      pixel[i] = (p1[i] - p0[i])*k;

  } else {

    double *p = imgPixel(areaImg, x1, y1);

    for(int i = 0; i < 4; ++i)

      pixel[i] = p[i]*k;

  }

}

void filterBW(Img *img, int size, double k) {

  if (!img->data) return;

  Img area = {};

  filterArea(&area, img);

  char *res = malloc(img->w * img->h);

  double sum[2][4] = {};

  int count[2] = {};

  k *= 3;

  double avg[4] = {};

  char *r = res;

  for(int y = 0; y < img->h; ++y) {

    for(int x = 0; x < img->w; ++x) {

      pixelAreaAvg(avg, &area, x - size, y - size, x + size, y + size);

      double a = avg[0] + avg[1] + avg[2];

      double *c = imgPixel(img, x, y);

      double cb = c[0] + c[1] + c[2];

      int value = (cb > a + k);

      for(int i = 0; i < 4; ++i)

        sum[value][i] += c[i];

      ++count[value];

      *r++ = value;

    }

  }

  imgDestroy(&area);

  for(int j = 0; j < 2; ++j)

    if (count[j])

      for(int i = 0; i < 4; ++i)

        sum[j][i] /= count[j];

  double *p = img->data;

  for(char *r = res, *e = r + img->w*img->h; r < e; ++r, p += 4)

    memcpy(p, sum[(int)*r], sizeof(*sum));

  free(res);

}

int bwGenerate(int argc, char **argv) {

  char *inFile = "", *outFile = "";

  double sub = 1.0/32;

  double k = -0.1;

  ArgDesc descs[] = {

    { 'i', &inFile, &pargString, "path to input image" },

    { 'o', &outFile, &pargString, "path to output image" },

    { 's', &sub, &pargDouble, "subrect size (relative to image size, default: 0.03)" },

    { 'k', &k, &pargDouble, "threshold shift (default: -0.1)" },

    {}

  };

  if (!parseArgs(descs, argc, argv) || !inFile[0] || !outFile[0]) {

    printUsage(descs);

    return FALSE;

  }

  Img img = {};

  if (!imgLoad(&img, inFile))

    return FALSE;

  filterBW(&img, (img.w + img.h)*sub/4, k);

  if (!imgSave(&img, outFile)) {

    imgDestroy(&img);

    return FALSE;

  }

  imgDestroy(&img);

  return TRUE;

}