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;
}