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

  }

}