double calcCurvePointEx(double *curveX, int stepX, double *curveY, int stepY, int count, double x) {
  if (count <= 0) return 0;
  if (count == 1) return *curveY;
  if (x <= curveX[0]) return *curveY;
  if (x >= curveX[(count-1)*stepX]) return curveY[(count-1)*stepY];

  int i = 0;
  while(i+1 < count && curveX[(i+1)*stepX] < x)
    ++i;

  int i0 = i > 0 ? i-1 : i;
  int i1 = i;
  int i2 = i1+1 < count ? i1+1 : i1;
  int i3 = i2+1 < count ? i2+1 : i2;

  double p0[] = { curveX[i0*stepX], curveY[i0*stepY] };
  double p1[] = { curveX[i1*stepX], curveY[i1*stepY] };
  double p2[] = { curveX[i2*stepX], curveY[i2*stepY] };
  double p3[] = { curveX[i3*stepX], curveY[i3*stepY] };

  double dx = p2[0] - p1[0];
  if (dx <= IMG_PRECISION) return p1[1];

  double a = p1[1];
  double b = p2[1];
  double ta = p2[0]-p0[0]; ta = ta > IMG_PRECISION ? (p2[1]-p0[1])/ta*dx : 0;
  double tb = p3[0]-p1[0]; tb = tb > IMG_PRECISION ? (p3[1]-p1[1])/tb*dx : 0;
  double l = (x - p1[0])/dx;

  return a + (ta + (3*(b-a)-ta-ta-tb + (2*(a-b)+ta+tb)*l)*l)*l;
}


double calcCurvePoint(double *curve, int count, double x)
  { return calcCurvePointEx(curve, 2, curve+1, 2, count, x); }


double calcTCurvePointEx(
  double *curveX, int stepX,
  double *curveY, int stepY,
  double *curveT, int stepT,
  int count, double x)
{
  if (count <= 0) return 0;
  if (count == 1) return *curveY;
  if (x <= curveX[0]) return *curveY;
  if (x >= curveX[(count-1)*stepX]) return curveY[(count-1)*stepY];

  int i0 = 0, i1 = 1;
  while(i1+1 < count && curveX[i1*stepX] < x)
    ++i0, ++i1;

  double x0 = curveX[i0*stepX];
  double x1 = curveX[i1*stepX];
  double dx = x1 - x0;
  if (dx <= IMG_PRECISION) return curveY[i0*stepY];

  double a  = curveY[i0*stepY];
  double b  = curveY[i1*stepY];
  double ta = curveT[i0*stepT]*dx;
  double tb = curveT[i1*stepT]*dx;
  double l = (x - x0)/dx;

  return a + (ta + (3*(b-a)-ta-ta-tb + (2*(a-b)+ta+tb)*l)*l)*l;
}


double calcTCurvePoint(double *curve, int count, double x)
  { return calcTCurvePointEx(curve, 3, curve+1, 3, curve+2, 3, count, x); }



void filterRasterizeCurve(Img *img, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double v = calcCurvePoint(curve, count, l0 + (l1 - l0)*x/(img->w - 1));
    for(int y = 0; y < img->h; ++y) {
      double *p = imgPixel(img, x, y);
      p[0] = p[1] = p[2] = p[3] = v;
    }
  }
}


void filterRasterizeCurveCh(Img *img, int ch, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double v = calcCurvePoint(curve, count, l0 + (l1 - l0)*x/(img->w - 1));
    for(int y = 0; y < img->h; ++y)
      imgPixel(img, x, y)[ch] = v;
  }
}


void filterRasterizeCurveColor(Img *img, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double l = l0 + (l1 - l0)*x/(img->w - 1);
    double c[] = {
      calcCurvePointEx(curve, 5, curve+1, 5, count, l),
      calcCurvePointEx(curve, 5, curve+2, 5, count, l),
      calcCurvePointEx(curve, 5, curve+3, 5, count, l),
      calcCurvePointEx(curve, 5, curve+4, 5, count, l) };
    for(int y = 0; y < img->h; ++y)
      memcpy(imgPixel(img, x, y), c, sizeof(c));
  }
}


void filterRasterizeTCurve(Img *img, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double v = calcTCurvePoint(curve, count, l0 + (l1 - l0)*x/(img->w - 1));
    for(int y = 0; y < img->h; ++y) {
      double *p = imgPixel(img, x, y);
      p[0] = p[1] = p[2] = p[3] = v;
    }
  }
}


void filterRasterizeTCurveCh(Img *img, int ch, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double v = calcTCurvePoint(curve, count, l0 + (l1 - l0)*x/(img->w - 1));
    for(int y = 0; y < img->h; ++y)
      imgPixel(img, x, y)[ch] = v;
  }
}


void filterRasterizeTCurveColor(Img *img, double *curve, int count, double l0, double l1) {
  if (!img->data) return;
  for(int x = 0; x < img->w; ++x) {
    double l = l0 + (l1 - l0)*x/(img->w - 1);
    double c[] = {
      calcTCurvePointEx(curve, 6, curve+1, 6, curve+2, 6, count, l),
      calcTCurvePointEx(curve, 6, curve+1, 6, curve+3, 6, count, l),
      calcTCurvePointEx(curve, 6, curve+1, 6, curve+4, 6, count, l),
      calcTCurvePointEx(curve, 6, curve+1, 6, curve+5, 6, count, l) };
    for(int y = 0; y < img->h; ++y)
      memcpy(imgPixel(img, x, y), c, sizeof(c));
  }
}