diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..f7b181e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+/data/output/
+/icy
+/icy-png
+/icy-heli
diff --git a/arg.inc.c b/arg.inc.c
new file mode 100644
index 0000000..b52ff5f
--- /dev/null
+++ b/arg.inc.c
@@ -0,0 +1,56 @@
+
+struct ArgDesc;
+typedef int (*ArgFunc)(const struct ArgDesc *desc, const char *arg);
+
+typedef struct ArgDesc {
+  char c;
+  void *data;
+  ArgFunc func;
+  const char *desc;
+} ArgDesc;
+
+
+int pargInt(const ArgDesc *desc, const char *arg)
+  { *(int*)desc->data = atof(arg); return TRUE; }
+int pargDouble(const ArgDesc *desc, const char *arg)
+  { *(double*)desc->data = atoi(arg); return TRUE; }
+int pargString(const ArgDesc *desc, const char *arg)
+  { *(const char**)desc->data = arg; return TRUE; }
+
+
+int parseArgs(const ArgDesc *descs, int argc, char **argv) {
+  for(; argc; --argc, ++argv) {
+    int found = FALSE;
+    for(const ArgDesc *desc = descs; desc->c && desc->data; ++desc) {
+      if ((*argv)[0] == '-' && (*argv)[1] == desc->c && !(*argv)[2]) {
+        if (desc->func) {
+          --argc; ++argv;
+          if (!argc || !desc->func(desc, *argv)) {
+            fprintf(stderr, "Wrong argument for -%c: %s\n", desc->c, *argv);
+            fflush(stderr);
+            return FALSE;
+          }
+        } else {
+          *((int*)desc->data) = TRUE;
+        }
+        found = TRUE;
+        break;
+      }
+    }
+    if (!found) {
+      fprintf(stderr, "Unknown arg: %s\n", *argv);
+      fflush(stderr);
+      return FALSE;
+    }
+  }
+  return TRUE;
+}
+
+
+void printUsage(const ArgDesc *descs) {
+  printf("Available command-line arguments:\n");
+  for(const ArgDesc *desc = descs; desc->c && desc->data; ++desc)
+    printf("   -%c %s - %s\n", desc->c, desc->func ? "<arg>" : "     ", desc->desc);
+  fflush(stdout);
+}
+
diff --git a/build-icy-heli.sh b/build-icy-heli.sh
new file mode 100755
index 0000000..047d25c
--- /dev/null
+++ b/build-icy-heli.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+
+set -e
+
+cc -Wall -DNDEBUG -O3 icy-heli.c -lm `pkg-config --cflags --libs helianthus` -o icy-heli
+
diff --git a/build-icy-png.sh b/build-icy-png.sh
new file mode 100755
index 0000000..4043ace
--- /dev/null
+++ b/build-icy-png.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+
+set -e
+
+cc -Wall -DNDEBUG -O3 icy-png.c -lm `pkg-config --cflags --libs libpng` -o icy-png
+
diff --git a/build-icy.sh b/build-icy.sh
new file mode 100755
index 0000000..b67e0ef
--- /dev/null
+++ b/build-icy.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+
+set -e
+
+cc -Wall -DNDEBUG -O3 icy.c -lm -o icy
+
+
+
diff --git a/data/icy-color.png b/data/icy-color.png
new file mode 100644
index 0000000..f7c260c
Binary files /dev/null and b/data/icy-color.png differ
diff --git a/data/icy-color.tga b/data/icy-color.tga
new file mode 100644
index 0000000..4c3bf08
Binary files /dev/null and b/data/icy-color.tga differ
diff --git a/data/icy-curve.png b/data/icy-curve.png
new file mode 100644
index 0000000..f995a02
Binary files /dev/null and b/data/icy-curve.png differ
diff --git a/data/icy-curve.tga b/data/icy-curve.tga
new file mode 100644
index 0000000..1300b9e
Binary files /dev/null and b/data/icy-curve.tga differ
diff --git a/data/icy-demo.png b/data/icy-demo.png
new file mode 100644
index 0000000..82765ad
Binary files /dev/null and b/data/icy-demo.png differ
diff --git a/data/input/text-large.png b/data/input/text-large.png
new file mode 100644
index 0000000..ab26642
Binary files /dev/null and b/data/input/text-large.png differ
diff --git a/data/input/text-large.tga b/data/input/text-large.tga
new file mode 100644
index 0000000..c4faddb
Binary files /dev/null and b/data/input/text-large.tga differ
diff --git a/data/input/text.png b/data/input/text.png
new file mode 100644
index 0000000..c9eb7bf
Binary files /dev/null and b/data/input/text.png differ
diff --git a/data/input/text.tga b/data/input/text.tga
new file mode 100644
index 0000000..ab12a62
Binary files /dev/null and b/data/input/text.tga differ
diff --git a/data/output/.placeholder b/data/output/.placeholder
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/data/output/.placeholder
diff --git a/filter.base.inc.c b/filter.base.inc.c
new file mode 100644
index 0000000..5ef0dc7
--- /dev/null
+++ b/filter.base.inc.c
@@ -0,0 +1,134 @@
+
+
+
+void filterFill(Img *img, double *pix) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4)
+    memcpy(p, pix, sizeof(double)*4);
+}
+
+
+void filterFillCh(Img *img, int ch, double val) {
+  if (!img->data) return;
+  for(double *p = img->data + ch, *e = p + 4*img->w*img->h; p < e; p += 4)
+    *p = val;
+}
+
+
+void filterMin(Img *img, double *pix) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    if (!(p[0] >= pix[0])) p[0] = pix[0];
+    if (!(p[1] >= pix[1])) p[1] = pix[1];
+    if (!(p[2] >= pix[2])) p[2] = pix[2];
+    if (!(p[3] >= pix[3])) p[3] = pix[3];
+  }
+}
+void filterMinCh(Img *img, int ch, double val) {
+  if (!img->data) return;
+  for(double *p = img->data + ch, *e = p + 4*img->w*img->h; p < e; p += 4)
+    if (!(*p <= val)) *p = val;
+}
+
+
+void filterMax(Img *img, double *pix) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    if (!(p[0] <= pix[0])) p[0] = pix[0];
+    if (!(p[1] <= pix[1])) p[1] = pix[1];
+    if (!(p[2] <= pix[2])) p[2] = pix[2];
+    if (!(p[3] <= pix[3])) p[3] = pix[3];
+  }
+}
+void filterMaxCh(Img *img, int ch, double val) {
+  if (!img->data) return;
+  for(double *p = img->data + ch, *e = p + 4*img->w*img->h; p < e; p += 4)
+    if (!(*p <= val)) *p = val;
+}
+
+
+void filterClampEx(Img *img, double *pmin, double *pmax) {
+  filterMax(img, pmax);
+  filterMin(img, pmin);
+}
+void filterClamp(Img *img) {
+  double pmin[4] = {}, pmax[4] = { 1, 1, 1, 1 };
+  filterClampEx(img, pmin, pmax);
+}
+void filterClampCh(Img *img, int ch, double vmin, double vmax) {
+  filterMaxCh(img, ch, vmax);
+  filterMinCh(img, ch, vmin);
+}
+
+
+void filterMulAlpha(Img *img) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    double a = p[3];
+    p[0] *= a; p[1] *= a; p[2] *= a;
+  }
+}
+
+
+void filterDivAlpha(Img *img) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    double k = p[3];
+    k = fabs(k) > IMG_PRECISION ? 1/k : 0;
+    p[0] *= k; p[1] *= k; p[2] *= k;
+  }
+}
+
+
+void filterResampleCubic(Img *img, int w, int h) {
+  Img imgC = {};
+  imgInit(&imgC, w, h);
+  imgSwap(img, &imgC);
+  if (img->data && imgC.data) {
+    double kx = (double)imgC.w/img->w;
+    double ky = (double)imgC.h/img->h;
+    for(int y = 0; y < img->h; ++y)
+      for(int x = 0; x < img->w; ++x)
+        imgCubic(&imgC, x*kx, y*ky, imgPixel(img, x, y));
+  }
+  imgDestroy(&imgC);
+}
+
+
+void filterCurve(Img *img, Img *icurve) {
+  if (!img->data || !icurve->data) return;
+  double k = icurve->w;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    p[0] = imgCubicHCh(icurve, 0, p[0]*k, 0);
+    p[1] = imgCubicHCh(icurve, 1, p[1]*k, 0);
+    p[2] = imgCubicHCh(icurve, 2, p[2]*k, 0);
+    p[3] = imgCubicHCh(icurve, 3, p[3]*k, 0);
+  }
+}
+
+
+void filterCurveCh(Img *img, Img *icurve, int dstCh, int srcCh, int curveCh) {
+  if (!img->data || !icurve->data) return;
+  double k = icurve->w;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4)
+    p[dstCh] = imgCubicHCh(icurve, curveCh, p[srcCh]*k, 0);
+}
+
+
+void filterCurveByCh(Img *img, Img *icurve, int ch) {
+  if (!img->data || !icurve->data) return;
+  double k = icurve->w;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4)
+    imgCubicH(icurve, p[ch]*k, 0, p);
+}
+
+
+void filterScreenToCompositeHSV(Img *img) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    rgbToHsv(p, p);
+    p[3] = p[2]; p[2] = 1;
+    hsvToRgb(p, p);
+  }
+}
+
diff --git a/filter.blend.inc.c b/filter.blend.inc.c
new file mode 100644
index 0000000..a59c47a
--- /dev/null
+++ b/filter.blend.inc.c
@@ -0,0 +1,106 @@
+
+
+typedef void (*BlendFunc)(double *dst, double *src, double val);
+
+
+void bfuncCopy(double *dst, double *src, double val)
+  { memcpy(dst, src, 4*sizeof(double)); }
+
+void bfuncMix(double *dst, double *src, double val) {
+  double iv = 1 - val;
+  dst[0] = dst[0]*iv + src[0]*val;
+  dst[1] = dst[1]*iv + src[1]*val;
+  dst[2] = dst[2]*iv + src[2]*val;
+  dst[3] = dst[3]*iv + src[3]*val;
+}
+
+void bfuncComposite(double *dst, double *src, double val) {
+  val *= src[3];
+  double iv = 1 - val;
+  dst[0] = dst[0]*iv + src[0]*val;
+  dst[1] = dst[1]*iv + src[1]*val;
+  dst[2] = dst[2]*iv + src[2]*val;
+  dst[3] = 1 - (1 - dst[3])*val;
+}
+
+void bfuncSum(double *dst, double *src, double val) {
+  val *= src[3];
+  dst[0] = dst[0] + src[0]*val;
+  dst[1] = dst[1] + src[1]*val;
+  dst[2] = dst[2] + src[2]*val;
+}
+
+void bfuncMul(double *dst, double *src, double val) {
+  val *= src[3];
+  double iv = 1 - val;
+  dst[0] = dst[0]*iv + dst[0]*src[0]*val;
+  dst[1] = dst[1]*iv + dst[1]*src[1]*val;
+  dst[2] = dst[2]*iv + dst[2]*src[2]*val;
+}
+
+
+
+void filterBlendEx(Img *img, Img *src, BlendFunc func, double val, int dx, int dy, int sx, int sy, int w, int h) {
+  if (!img->data || !src->data || !func) return;
+
+  if (sx < 0) { dx -= sx; w += sx; sx = 0; }
+  if (sy < 0) { dy -= sy; h += sy; sy = 0; }
+  if (dx < 0) { sx -= dx; w += dx; dx = 0; }
+  if (dy < 0) { sy -= dy; h += dy; dy = 0; }
+
+  if (w > img->w - dx) w = img->w - dx;
+  if (h > img->h - dy) h = img->h - dy;
+  if (w > src->w - sx) w = src->w - sx;
+  if (h > src->h - sy) h = src->h - sy;
+
+  if (w <= 0 || h <= 0) return;
+  int dr = 4*(img->w - w);
+  int sr = 4*(src->w - w);
+
+  if (img == src) {
+    int offset = ((sy - dy)*img->w + sx - dx)*4;
+    if (offset > 0) {
+      for(double *pd = imgPixel(img, dx, dy), *e = pd + 4*h*img->w; pd < e; pd += dr)
+        for(double *re = pd + 4*w; pd < re; pd += 4)
+          func(pd, pd + offset, val);
+    } else
+    if (offset < 0) {
+      for(double *pd = imgPixel(img, dx+w-1, dy+h-1), *e = pd - 4*h*img->w; pd > e; pd -= dr)
+        for(double *re = pd - 4*w; pd > re; pd -= 4)
+          func(pd, pd + offset, val);
+    } else {
+      double ps[4];
+      for(double *pd = imgPixel(img, dx, dy), *e = pd + 4*h*img->w; pd < e; pd += dr)
+        for(double *re = pd + 4*w; pd < re; pd += 4)
+          { memcpy(ps, pd, sizeof(ps)); func(pd, ps, val); }
+    }
+  } else {
+    for(double *pd = imgPixel(img, dx, dy), *ps = imgPixel(src, sx, sy), *e = pd + 4*h*img->w; pd < e; pd += dr, ps += sr)
+      for(double *re = pd + 4*w; pd < re; pd += 4, ps += 4)
+        func(pd, ps, val);
+  }
+}
+
+void filterBlend(Img *img, Img *src, BlendFunc func, double val)
+  { filterBlendEx(img, src, func, val, 0, 0, 0, 0, src->w, src->h); }
+
+
+
+
+void imgCopyTo(Img *img, Img *src, int dx, int dy, int sx, int sy, int w, int h)
+  { filterBlendEx(img, src, bfuncCopy, 1, dx, dy, sx, sy, w, h); }
+
+void imgCropTo(Img *img, Img *src, int x, int y, int w, int h) {
+  imgInit(img, w, h);
+  imgCopyTo(img, src, 0, 0, x, y, w, h);
+}
+
+void imgCrop(Img *img, int x, int y, int w, int h) {
+  Img imgC = {};
+  imgSwap(img, &imgC);
+  imgCropTo(img, &imgC, x, y, w, h);
+  imgDestroy(&imgC);
+}
+
+void imgCopy(Img *img, Img *src)
+  { imgCropTo(img, src, 0, 0, src->w, src->h); }
diff --git a/filter.blur.inc.c b/filter.blur.inc.c
new file mode 100644
index 0000000..d8c63bc
--- /dev/null
+++ b/filter.blur.inc.c
@@ -0,0 +1,196 @@
+
+
+int toPowerOf2(int x) {
+	int xx;
+  for(xx = 1; x > 1; x >>= 1)
+    xx <<= 1;
+  return xx;
+}
+
+
+
+void fft(double *real, int realStep, double *img, int imgStep, int count, int invert) {
+	if (count < 2 || !realStep || !imgStep) return;
+  count = toPowerOf2(count);
+
+  // check order
+  int flipReal = 0;
+  if (realStep < 0) {
+    flipReal = 1;
+    real += realStep*(count - 1);
+    realStep = -realStep;
+  }
+  int flipImg = 0;
+  if (imgStep < 0) {
+    flipImg = 1;
+    img += imgStep*(count - 1);
+    imgStep = -imgStep;
+  }
+
+	double *realEnd = real + count*realStep;
+	double *imgEnd  = img  + count*imgStep;
+
+  // flip if step was negative
+  if (flipReal)
+    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
+      { double x = *i; *i = *j; *j = x; }
+  if (flipImg)
+    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
+      { double x = *i; *i = *j; *j = x; }
+
+
+	// bit-reversal permutation
+	for(int i=0, j=0; i < count; ++i) {
+		if (j > i) {
+      double x, *a, *b;
+      a = real + i*realStep; b = real + j*realStep; x = *a; *a = *b; *b = x;
+      a = img  + i*imgStep;  b = img  + j*imgStep;  x = *a; *a = *b; *b = x;
+    }
+		int m = count/2;
+		while(m >= 1 && j >= m)
+			{ j -= m; m /= 2; }
+		j += m;
+	}
+
+	double k = invert ? PI : -PI;
+	for(int mmax = 1; mmax < count; mmax *= 2) {
+		// rotation coefficients
+		double angle = k/mmax;
+		double wpR = sin(0.5*angle), wpI = sin(angle);
+		wpR *= wpR*2;
+
+		double wR = 1, wI = 0;
+		int mmaxRealStep  = mmax*realStep;
+		int mmaxRealStep2 = 2*mmaxRealStep;
+		int mmaxImgStep   = mmax*imgStep;
+		int mmaxImgStep2  = 2*mmaxImgStep;
+		for(int m = 0; m < mmax; ++m) {
+			// rotate w
+			double wwR = wR;
+			wR = wR - wwR*wpR - wI*wpI,
+      wI = wI + wwR*wpI - wI*wpR;
+			// process subsequences
+			for( double *iR = real + m*realStep,
+			            *iI = img  + m*imgStep,
+			            *jR = iR + mmaxRealStep,
+                  *jI = iI + mmaxImgStep;
+           iR < realEnd;
+           iR += mmaxRealStep2,
+           iI += mmaxImgStep2,
+           jR += mmaxRealStep2,
+           jI += mmaxImgStep2 )
+      {
+				// radix
+				double tR = *jR*wR - *jI*wI;
+				double tI = *jR*wI + *jI*wR;
+				*jR = *iR - tR;
+				*jI = *iI - tI;
+				*iR += tR;
+				*iI += tI;
+			}
+		}
+	}
+
+	// reverse order
+  if (!flipReal)
+    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
+      { double x = *i; *i = *j; *j = x; }
+  if (!flipImg)
+    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
+      { double x = *i; *i = *j; *j = x; }
+
+	// divide by count to complete back-FFT
+	if (invert) {
+		double k = 1.0/count;
+		for(double *i = real; i < realEnd; i += realStep)
+			*i *= k;
+		for(double *i = img; i < imgEnd; i += imgStep)
+			*i *= k;
+	}
+}
+
+
+void fft2d(double *real, int realColStep, int realRowStep, double *img, int imgColStep, int imgRowStep, int rows, int cols, int invert) {
+	if (rows < 1 || cols < 1 || !realColStep || !imgColStep || !realRowStep || imgRowStep) return;
+  rows = toPowerOf2(rows);
+  cols = toPowerOf2(cols);
+
+	// fft rows
+	if (cols > 1)
+		for(double *r = real, *i = img, *end = real + rows*realRowStep; r < end; r += realRowStep, i += imgRowStep)
+			fft(r, realColStep, i, imgColStep, cols, invert);
+	// fft cols
+	if (rows > 1)
+		for(double *r = real, *i = img, *end = real + cols*realColStep; r < end; r += realColStep, i += imgColStep)
+			fft(r, realRowStep, i, imgRowStep, rows, invert);
+}
+
+
+double preGauss(double x, double r) {
+  if (fabs(r) < IMG_PRECISION)
+    return fabs(x) < IMG_PRECISION ? 1 : 0;
+  return exp(-0.5*x*x/(r*r));
+}
+
+
+double gauss(double x, double r) {
+  static const double k = 1/sqrt(2*PI);
+  return k*preGauss(x, r)/fabs(r);
+}
+
+
+void gaussBlurRows(double *row, int colStep, int rowStep, int cols, int rows, double size) {
+  if (!(size > IMG_PRECISION)) return;
+  int cnt = cols + 2*(int)ceil(size*4 - IMG_PRECISION);
+  int p2cnt = toPowerOf2(cnt);
+  if (p2cnt < cnt) p2cnt <<= 1;
+  cnt = p2cnt;
+
+  double *buf = (double*)calloc(1, sizeof(double)*4*cnt);
+  double *pat = buf + 2*cnt;
+  pat[0]   = preGauss(0,     size);
+  pat[cnt] = preGauss(cnt/2, size);
+  for(int i = 2; i < cnt; i += 2)
+    pat[i] = pat[cnt*2 - i] = preGauss(i, size);
+  fft(pat, 2, pat+1, 2, cnt, FALSE);
+
+  double *buf0 = buf + (cnt - cols);
+  double *buf1 = buf0 + (cols - 1)*2;
+  for(int i = 0; i < rows; ++i) {
+    for(int ch = 0; ch < 4; ++ch) {
+      memset(buf, 0, sizeof(double)*2*cnt);
+      double *chrow = row + i*rowStep + ch;
+      for(int j = 0; j < cols; ++j)
+        buf0[j*2] = chrow[j*colStep];
+      for(double v = *buf0, *b = buf;      b < buf0; b += 2) *b = v;
+      for(double v = *buf1, *b = buf1 + 2; b < pat;  b += 2) *b = v;
+      fft(buf, 2, buf+1, 2, cnt, FALSE);
+      for(double *b = buf, *p = pat; b < pat; b += 2, p += 2) {
+        double bR = b[0], bI = b[1];
+        b[0] = bR*p[0] - bI*p[1];
+        b[1] = bR*p[1] + bI*p[0];
+      }
+      fft(buf, 2, buf+1, 2, cnt, TRUE);
+      for(int j = 0; j < cols; ++j) {
+        double *b = buf0 + j*2;
+        chrow[j*colStep] = sqrt(b[0]*b[0] + b[1]*b[1]);
+      }
+    }
+  }
+  free(buf);
+}
+
+
+void filterBlurGauss(Img *img, double sx, double sy) {
+  if (!img->data) return;
+  sx = fabs(sx);
+  sy = fabs(sy);
+  if (!(sx > IMG_PRECISION) && !(sy > IMG_PRECISION)) return;
+
+  filterMulAlpha(img);
+  gaussBlurRows(img->data, 4, 4*img->w, img->w, img->h, sx);
+  gaussBlurRows(img->data, 4*img->w, 4, img->h, img->w, sy);
+  filterDivAlpha(img);
+}
+
+
diff --git a/filter.icy.inc.c b/filter.icy.inc.c
new file mode 100644
index 0000000..f29fb62
--- /dev/null
+++ b/filter.icy.inc.c
@@ -0,0 +1,133 @@
+
+
+void filterIcy(Img *img, Img *icyCurve, Img *icyColor, double size, int withWind) {
+  if (!img->data) return;
+  size = fabs(size);
+  if (!(size > IMG_PRECISION)) withWind = FALSE;
+
+  const int    optNoiseWidth0 = (int)round(size/4);
+  const int    optNoiseWidth1 = (int)round(size);
+  const double optNoise0 = 1;
+  const double optNoise1 = 0.75;
+  const double optBright0 = 1.1;
+  const double optBright1 = 1.4;
+
+  const double optWind = 0.75;
+  const int    optWindWidth0 = (int)round(size/4);
+  const int    optWindWidth1 = (int)round(size/2);
+  const double optWindNoise0 = 0.5;
+  const double optWindNoise1 = 0.75;
+  const double optWindAtt = withWind ? pow(0.95, 4/size) : 0;
+
+  double blurSize = size/2;
+  if (optNoiseWidth0) filterMulNoiseGray(img, img->w/optNoiseWidth0, img->h/optNoiseWidth0, 1-optNoise0, 1);
+  if (optNoiseWidth1) filterMulNoiseGray(img, img->w/optNoiseWidth1, img->h/optNoiseWidth1, 1-optNoise1, 1);
+  filterFillCh(img, 3, 1);
+  filterBlurGauss(img, blurSize, blurSize);
+
+  ClMat m;
+  clmMulColor(&m, optBright0);
+  filterMatrix(img, &m);
+
+  if (icyCurve && icyCurve->data) {
+    //int mi = 0;
+    //for(int i = 1; i < icyCurve->w; ++i)
+    //  if (imgPixel(icyCurve, i-1, 0)[0] - IMG_PRECISION <= imgPixel(icyCurve, i, 0)[0])
+    //    mi = i; else break;
+    int mi = 128;
+    int mmi = size > 4 ? round(mi*4/size) : mi;
+    if (mmi < mi && mi < icyCurve->w - 1) {
+      Img imgA = {}, imgB = {};
+      imgCropTo(&imgA, icyCurve, 0, 0, mi, icyCurve->h);
+      filterResampleCubic(&imgA, mmi, icyCurve->h);
+      imgInit(&imgB, icyCurve->w - mi + mmi + 1, icyCurve->h);
+      imgCopyTo(&imgB, &imgA, 1, 0, 0, 0, imgA.w, imgA.h);
+      imgCopyTo(&imgB, icyCurve, mmi + 1, 0, mi, 0, icyCurve->w - mi, icyCurve->h);
+      filterFillCh(&imgB, 3, 1);
+      filterCurveByCh(img, &imgB, 0);
+      imgDestroy(&imgA);
+      imgDestroy(&imgB);
+    } else {
+      filterCurveByCh(img, icyCurve, 0);
+    }
+  }
+
+  clmMulColor(&m, optBright1);
+  filterMatrix(img, &m);
+
+  if (withWind) {
+    Img wind = {};
+    imgCopy(&wind, img);
+
+    Img pat = {};
+    imgInit(&pat, img->w, 1);
+    double clOne[] = {1,1,1,1};
+    filterFill(&pat, clOne);
+    if (optWindWidth0) filterMulNoiseGray(&pat, img->w/optWindWidth0, 1, 1-optWindNoise0, 1);
+    if (optWindWidth1) filterMulNoiseGray(&pat, img->w/optWindWidth1, 1, 1-optWindNoise1, 1);
+    filterWind(&wind, &pat, optWindAtt);
+    imgDestroy(&pat);
+
+    filterBlend(&wind, img, &bfuncSum, -1);
+    filterBlend(img, &wind, &bfuncSum, optWind);
+    imgDestroy(&wind);
+  }
+
+  if (icyColor)
+    filterCurveByCh(img, icyColor, 0);
+
+  //clmScreenToCompositeYUV(&m);
+  //filterMatrix(img, &m);
+  filterScreenToCompositeHSV(img);
+}
+
+
+int icyGenerate(Img *img, int argc, char **argv) {
+  if (img) imgDestroy(img);
+
+  char *inFile = "", *outFile = "";
+  char *curveFile = "data/icy-curve.tga", *colorFile = "data/icy-color.tga";
+  double size = 4;
+  int withWind = TRUE;
+  int withoutWind = FALSE;
+  ArgDesc descs[] = {
+    { 'i', &inFile, &pargString, "path to input image" },
+    { 'o', &outFile, &pargString, "path to output image" },
+    { 'c', &curveFile, &pargString, "path to image with crystal curve" },
+    { 'l', &colorFile, &pargString, "path to image with color curve" },
+    { 's', &size, &pargDouble, "effect size" },
+    { 'w', &withWind, NULL, "add icicles (default)" },
+    { 'n', &withoutWind, NULL, "no icicles" },
+    {}
+  };
+
+  if (!parseArgs(descs, argc, argv) || !inFile[0] || !outFile[0]) {
+    printUsage(descs);
+    return FALSE;
+  }
+
+  if (withoutWind) withWind = FALSE;
+
+  Img imgC = {};
+  Img imgCurve = {};
+  Img imgColor = {};
+  if (!imgLoad(&imgC, inFile))
+    return FALSE;
+
+  imgLoad(&imgCurve, curveFile);
+  imgLoad(&imgColor, colorFile);
+  filterIcy(&imgC, &imgCurve, &imgColor, size, withWind);
+
+  if (!imgSave(&imgC, outFile))
+    imgDestroy(&imgC);
+  int success = !!imgC.data;
+
+  if (img) imgSwap(img, &imgC);
+  imgDestroy(&imgCurve);
+  imgDestroy(&imgCurve);
+  imgDestroy(&imgC);
+
+  return success;
+}
+
+
diff --git a/filter.matrix.inc.c b/filter.matrix.inc.c
new file mode 100644
index 0000000..bcd5d75
--- /dev/null
+++ b/filter.matrix.inc.c
@@ -0,0 +1,119 @@
+
+
+
+typedef struct {
+  double m[5][4];
+} ClMat;
+
+typedef void (*ClmInitFunc)(ClMat *m);
+
+
+void clmZero(ClMat *m) {
+  for(int r = 0; r < 5; ++r)
+    for(int c = 0; c < 4; ++c)
+      m->m[r][c] = 0;
+}
+
+
+void clmIdentity(ClMat *m) {
+  for(int r = 0; r < 5; ++r)
+    for(int c = 0; c < 4; ++c)
+      m->m[r][c] = r == c;
+}
+
+
+void clmCopy(ClMat *dst, ClMat *src)
+  { if (dst != src) memcpy(dst, src, sizeof(*dst)); }
+
+
+void clmMulColor(ClMat *m, double val) {
+  for(int r = 0; r < 5; ++r)
+    for(int c = 0; c < 4; ++c)
+      m->m[r][c] = r == c ? (r < 3 ? val : 1) : 0;
+}
+
+
+void clmMulPix(double *dst, const ClMat *m, const double *pix) {
+  if (dst == pix) {
+    double p[4];
+    memcpy(p, pix, sizeof(p));
+    clmMulPix(dst, m, p);
+    return;
+  }
+
+  dst[0] = m->m[0][0]*pix[0] + m->m[1][0]*pix[1] + m->m[2][0]*pix[2] + m->m[3][0]*pix[3] + m->m[4][0];
+  dst[1] = m->m[0][1]*pix[0] + m->m[1][1]*pix[1] + m->m[2][1]*pix[2] + m->m[3][1]*pix[3] + m->m[4][1];
+  dst[2] = m->m[0][2]*pix[0] + m->m[1][2]*pix[1] + m->m[2][2]*pix[2] + m->m[3][2]*pix[3] + m->m[4][2];
+  dst[3] = m->m[0][3]*pix[0] + m->m[1][3]*pix[1] + m->m[2][3]*pix[2] + m->m[3][3]*pix[3] + m->m[4][3];
+}
+
+
+void clmMul(ClMat *dst, const ClMat *a, const ClMat *b) {
+  if (dst == a || dst == b) {
+    ClMat m;
+    clmCopy(&m, dst);
+    clmMul(dst, dst == a ? &m : a, dst == b ? &m : b);
+    return;
+  }
+
+  clmMulPix(dst->m[0], a, b->m[0]);
+  clmMulPix(dst->m[1], a, b->m[1]);
+  clmMulPix(dst->m[2], a, b->m[2]);
+  clmMulPix(dst->m[3], a, b->m[3]);
+  clmMulPix(dst->m[4], a, b->m[4]);
+}
+
+
+void clmGray(ClMat *m) {
+  for(int r = 0; r < 5; ++r)
+    for(int c = 0; c < 4; ++c)
+      m->m[r][c] = r < 3 && c < 3 ? 1/3.0 : r == c;
+}
+
+
+void clmToYUV(ClMat *m) {
+  m->m[0][0] = 0.299; m->m[0][1] = -0.168736; m->m[0][2] =  0.5;      m->m[0][3] = 0;
+  m->m[1][0] = 0.587; m->m[1][1] =  0.331264; m->m[1][2] = -0.418688; m->m[1][3] = 0;
+  m->m[2][0] = 0.114; m->m[2][1] =  0.5;      m->m[2][2] = -0.081312; m->m[2][3] = 0;
+  m->m[3][0] = 0;     m->m[3][1] =  0;        m->m[3][2] =  0;        m->m[3][3] = 1;
+  m->m[4][0] = 0;     m->m[4][1] =  0;        m->m[4][2] =  0;        m->m[4][3] = 0;
+}
+
+
+void clmFromYUV(ClMat *m) {
+  m->m[0][0] = 1;     m->m[0][1] =  1;        m->m[0][2] = 1;     m->m[0][3] = 0;
+  m->m[1][0] = 0;     m->m[1][1] = -0.344136; m->m[1][2] = 1.772; m->m[1][3] = 0;
+  m->m[2][0] = 0.402; m->m[2][1] = -0.714136; m->m[2][2] = 0;     m->m[2][3] = 0;
+  m->m[3][0] = 0;     m->m[3][1] =  0;        m->m[3][2] = 0;     m->m[3][3] = 1;
+  m->m[4][0] = 0;     m->m[4][1] =  0;        m->m[4][2] = 0;     m->m[4][3] = 0;
+}
+
+
+void clmScreenToCompositeYUV(ClMat *m) {
+  // to YUV
+  clmToYUV(m);
+  // Alpha = Y-Luminance, Y-Luminance = 1
+  m->m[0][3] = m->m[0][0];
+  m->m[1][3] = m->m[1][0];
+  m->m[2][3] = m->m[2][0];
+  m->m[3][3] = m->m[3][0];
+  m->m[0][0] = 0;
+  m->m[1][0] = 0;
+  m->m[2][0] = 0;
+  m->m[4][0] = 1;
+  // YUV -> RGB
+  ClMat mm;
+  clmFromYUV(&mm);
+  clmMul(m, &mm, m);
+}
+
+
+
+
+void filterMatrix(Img *img, const ClMat *m) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4)
+    clmMulPix(p, m, p);
+}
+
+
diff --git a/filter.noise.inc.c b/filter.noise.inc.c
new file mode 100644
index 0000000..8b8821f
--- /dev/null
+++ b/filter.noise.inc.c
@@ -0,0 +1,52 @@
+
+
+void filterNoise(Img *img) {
+  if (!img->data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; p += 4) {
+    p[0] = randomFloat();
+    p[1] = randomFloat();
+    p[2] = randomFloat();
+    p[3] = randomFloat();
+  }
+}
+
+
+void filterMulNoise(Img *img, int w, int h, double vmin, double vmax) {
+  if (!img->data || w <= 0 || h <= 0) return;
+  Img imgC = {};
+  imgInit(&imgC, w, h);
+  filterNoise(&imgC);
+  for(int y = 0; y < img->h; ++y) {
+    for(int x = 0; x < img->w; ++x) {
+      double pn[4];
+      imgCubic(&imgC, (double)x/img->w*w, (double)y/img->h*h, pn);
+      double *p = imgPixel(img, x, y);
+      p[0] *= vmin + (vmax - vmin)*pn[0];
+      p[1] *= vmin + (vmax - vmin)*pn[1];
+      p[2] *= vmin + (vmax - vmin)*pn[2];
+      p[3] *= vmin + (vmax - vmin)*pn[3];
+    }
+  }
+  imgDestroy(&imgC);
+}
+
+
+void filterMulNoiseGray(Img *img, int w, int h, double vmin, double vmax) {
+  if (!img->data || w <= 0 || h <= 0) return;
+  Img imgC = {};
+  imgInit(&imgC, w, h);
+  filterNoise(&imgC);
+  for(int y = 0; y < img->h; ++y) {
+    for(int x = 0; x < img->w; ++x) {
+      double pn = vmin + (vmax - vmin)*imgCubicCh(&imgC, 0, (double)x/img->w*w, (double)y/img->h*h);
+      double *p = imgPixel(img, x, y);
+      p[0] *= pn;
+      p[1] *= pn;
+      p[2] *= pn;
+      //p[3] *= pn;
+    }
+  }
+  imgDestroy(&imgC);
+}
+
+
diff --git a/filter.wind.inc.c b/filter.wind.inc.c
new file mode 100644
index 0000000..6db58b7
--- /dev/null
+++ b/filter.wind.inc.c
@@ -0,0 +1,24 @@
+
+
+
+void filterWind(Img *img, Img *pattern, double att) {
+  if (!img->data) return;
+
+  for(int x = 0; x < img->w; ++x) {
+    double nk = pattern && pattern->data ? imgCubicHCh(pattern, 0, (double)x/img->w*pattern->w, 0) : 1;
+    double acc[4] = {};
+    double a = 1 - pow(1 - att, nk);
+    for(int y = 0; y < img->h; ++y) {
+      double *cl = imgPixel(img, x, y);
+      for(int c = 0; c < 4; ++c) {
+        if (cl[c] < acc[c]) {
+          cl[c] = acc[c];
+        } else
+        if (acc[c] < cl[c]) {
+          acc[c] = cl[c];
+        }
+        acc[c] *= a;
+      }
+    }
+  }
+}
diff --git a/icy-heli.c b/icy-heli.c
new file mode 100644
index 0000000..9758670
--- /dev/null
+++ b/icy-heli.c
@@ -0,0 +1,74 @@
+
+#include <math.h>
+#include <ctype.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <helianthus.h>
+
+
+#include "arg.inc.c"
+#include "img.inc.c"
+#include "img.ldr.inc.c"
+#include "filter.base.inc.c"
+#include "filter.matrix.inc.c"
+#include "filter.blend.inc.c"
+#include "filter.noise.inc.c"
+#include "filter.blur.inc.c"
+#include "filter.wind.inc.c"
+#include "filter.icy.inc.c"
+
+
+
+Animation anim;
+char *argvv[] = {"-i", "data/input/text.png", "-o", "data/output/icy-text.png", "-s", "4" };
+char **argv;
+int argc;
+
+
+void generate() {
+  if (anim) {
+    animationDestroy(anim);
+    anim = NULL;
+  }
+  Img img = {};
+  if (icyGenerate(&img, argc, argv)) {
+    unsigned char *data = imgToInt(&img);
+    anim = createAnimationFromImage(img.w, img.h, data, FALSE);
+    free(data);
+    if (img.w < 512) {
+      windowSetWidth(img.w);
+      windowSetHeight(img.h);
+    }
+    imgDestroy(&img);
+  }
+}
+
+
+void init() {
+  background(COLOR_BLACK);
+  generate();
+}
+
+
+void draw() {
+  if (keyWentDown("space")) generate();
+  noStroke();
+  rectTextured(anim, 0, 0, windowGetWidth(), windowGetHeight());
+}
+
+
+int main(int largc, char **largv) {
+  if (largc > 1) {
+    argc = largc - 1;
+    argv = largv + 1;
+  } else {
+    argc = sizeof(argvv)/sizeof(*argvv);
+    argv = argvv;
+  }
+  windowSetInit(&init);
+  windowSetDraw(&draw);
+  windowRun();
+  return 0;
+}
diff --git a/icy-png.c b/icy-png.c
new file mode 100644
index 0000000..701c6c4
--- /dev/null
+++ b/icy-png.c
@@ -0,0 +1,4 @@
+
+#include <png.h>
+#include "icy.c"
+
diff --git a/icy.c b/icy.c
new file mode 100644
index 0000000..c21be20
--- /dev/null
+++ b/icy.c
@@ -0,0 +1,39 @@
+
+#include <math.h>
+#include <ctype.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#ifndef TRUE
+#define TRUE 1
+#endif
+
+#ifndef FALSE
+#define FALSE 0
+#endif
+
+#ifndef PI
+#define PI 3.14159265358979323846
+#endif
+
+
+double randomFloat()
+ { return (double)rand()/RAND_MAX; }
+
+
+#include "arg.inc.c"
+#include "img.inc.c"
+#include "img.ldr.inc.c"
+#include "filter.base.inc.c"
+#include "filter.matrix.inc.c"
+#include "filter.blend.inc.c"
+#include "filter.noise.inc.c"
+#include "filter.blur.inc.c"
+#include "filter.wind.inc.c"
+#include "filter.icy.inc.c"
+
+
+int main(int argc, char **argv)
+  { return !icyGenerate(NULL, argc - 1, argv + 1); }
+
diff --git a/img.inc.c b/img.inc.c
new file mode 100644
index 0000000..e6421c3
--- /dev/null
+++ b/img.inc.c
@@ -0,0 +1,255 @@
+
+
+#define IMG_PRECISION (1e-8)
+
+
+typedef struct {
+  int w, h;
+  double *data;
+} Img;
+
+
+void setRgb(double *rgb, double r, double g, double b)
+  { rgb[0] = r; rgb[1] = g; rgb[2] = b; }
+
+int colorMin(const double *rgb)
+  { return rgb[0] < rgb[1] ? (rgb[0] < rgb[2] ? 0 : 2) : (rgb[1] < rgb[2] ? 1 : 2); }
+int colorMax(const double *rgb)
+  { return rgb[0] < rgb[1] ? (rgb[1] < rgb[2] ? 2 : 1) : (rgb[0] < rgb[2] ? 2 : 0); }
+
+void rgbToHsv(double *hsv, const double *rgb) {
+  int cmin = colorMin(rgb);
+  int cmax = colorMax(rgb);
+  double d = rgb[cmax] - rgb[cmin];
+
+  double h = 0;
+  if (d > IMG_PRECISION) {
+    double k = 1.0/d;
+    switch(cmax){
+    case 0: h = (rgb[1] - rgb[2])*k + 0; break;
+    case 1: h = (rgb[2] - rgb[0])*k + 2; break;
+    case 2: h = (rgb[0] - rgb[1])*k + 4; break;
+    }
+    h /= 6; h -= floor(h); h *= 360;
+  }
+  double s = rgb[cmax] > IMG_PRECISION ? d/rgb[cmax] : 0;
+  double v = rgb[cmax];
+  setRgb(hsv, h, s, v);
+}
+
+
+void hsvToRgb(double *rgb, const double *hsv) {
+  double h = hsv[0], s = hsv[1], v = hsv[2];
+  h -= floor(h/360)*360;
+  h /= 60.0;
+  int i = (int)h;
+  double f = h - i;
+  double p = v*(1 - s);
+  double q = v*(1 - s*f);
+  double t = v*(1 - s*(1 - f));
+  switch(i) {
+  case 0: setRgb(rgb, v, t, p); return;
+  case 1: setRgb(rgb, q, v, p); return;
+  case 2: setRgb(rgb, p, v, t); return;
+  case 3: setRgb(rgb, p, q, v); return;
+  case 4: setRgb(rgb, t, p, v); return;
+  case 5: setRgb(rgb, v, p, q); return;
+  }
+  return setRgb(rgb, v, t, p);
+}
+
+
+double cubicInterpolation(double p0, double p1, double p2, double p3, double l) {
+  double ll = l*l;
+  double lll = ll*l;
+  return 0.5*( p0*(  -lll + 2*ll - l)
+             + p1*( 3*lll - 5*ll + 2)
+             + p2*(-3*lll + 4*ll + l)
+             + p3*(   lll -   ll    ) );
+}
+
+double cubicRow(double *row, int count, double x) {
+  if (count < 0) return 0;
+  if (!(x > 0)) x = 0;
+  if (!(x < count)) x = count;
+  double xi = floor(x);
+  double l = x - xi;
+
+  int i = (int)x - 1;
+  double p0 = i < 0 ? row[0] : i >= count ? row[count-1] : row[i]; ++i;
+  double p1 = i < 0 ? row[0] : i >= count ? row[count-1] : row[i]; ++i;
+  double p2 = i < 0 ? row[0] : i >= count ? row[count-1] : row[i]; ++i;
+  double p3 = i < 0 ? row[0] : i >= count ? row[count-1] : row[i]; ++i;
+
+  return cubicInterpolation(p0, p1, p2, p3, l);
+}
+
+
+
+void imgDestroy(Img *img) {
+  free(img->data);
+  img->data = NULL;
+  img->w = img->h = 0;
+}
+
+
+void imgSwap(Img *imgA, Img *imgB) {
+  Img imgC;
+  memcpy(&imgC, imgA, sizeof(imgC));
+  memcpy(imgA, imgB, sizeof(imgC));
+  memcpy(imgB, &imgC, sizeof(imgC));
+}
+
+
+int imgInit(Img *img, int w, int h) {
+  imgDestroy(img);
+  if (w < 0 || h < 0) return FALSE;
+  img->data = calloc(1, sizeof(double)*4*w*h);
+  if (!img->data) return FALSE;
+  img->w = w;
+  img->h = h;
+  return TRUE;
+}
+
+
+double* imgPixel(const Img *img, int x, int y) {
+  if (!img->data) return NULL;
+  if (x >= img->w) x = img->w - 1;
+  if (y >= img->h) y = img->h - 1;
+  if (x < 0) x = 0;
+  if (y < 0) y = 0;
+  return img->data + (img->w*y + x)*4;
+}
+
+
+double* imgNearest(const Img *img, double x, double y)
+  { return imgPixel(img, (int)round(x), (int)round(y)); }
+
+
+double imgCubicHCh(const Img *img, int ch, double x, int y) {
+  if (!img->data) return 0;
+  if (!(x > 0)) x = 0;
+  if (!(x < img->w)) x = img->w;
+  double xi = floor(x);
+  int i = (int)xi - 1;
+  double p[] = {
+    imgPixel(img, i+0, y)[ch],
+    imgPixel(img, i+1, y)[ch],
+    imgPixel(img, i+2, y)[ch],
+    imgPixel(img, i+3, y)[ch] };
+  return cubicInterpolation(p[0], p[1], p[2], p[3], x - xi);
+}
+
+
+double imgCubicVCh(const Img *img, int ch, int x, double y) {
+  if (!img->data) return 0;
+  if (!(y > 0)) x = 0;
+  if (!(y < img->h)) y = img->h;
+  double yi = floor(y);
+  int i = (int)yi - 1;
+  double p[] = {
+    imgPixel(img, i+0, y)[ch],
+    imgPixel(img, i+1, y)[ch],
+    imgPixel(img, i+2, y)[ch],
+    imgPixel(img, i+3, y)[ch] };
+  return cubicInterpolation(p[0], p[1], p[2], p[3], y - yi);
+}
+
+
+void imgCubicH(const Img *img, double x, int y, double *pix) {
+  if (!img->data) { pix[0] = pix[1] = pix[2] = pix[3] = 0; return; }
+  if (!(x > 0)) x = 0;
+  if (!(x < img->w)) x = img->w;
+  double xi = floor(x);
+  double l = x - xi;
+
+  int i = (int)xi - 1;
+  double *p[] = {
+    imgPixel(img, i+0, y),
+    imgPixel(img, i+1, y),
+    imgPixel(img, i+2, y),
+    imgPixel(img, i+3, y) };
+
+  pix[0] = cubicInterpolation(p[0][0], p[1][0], p[2][0], p[3][0], l);
+  pix[1] = cubicInterpolation(p[0][1], p[1][1], p[2][1], p[3][1], l);
+  pix[2] = cubicInterpolation(p[0][2], p[1][2], p[2][2], p[3][2], l);
+  pix[3] = cubicInterpolation(p[0][3], p[1][3], p[2][3], p[3][3], l);
+}
+
+
+void imgCubicV(const Img *img, int x, double y, double *pix) {
+  if (!img->data) { pix[0] = pix[1] = pix[2] = pix[3] = 0; return; }
+  if (!(y > 0)) y = 0;
+  if (!(y < img->h)) y = img->h;
+  double yi = floor(y);
+  double l = y - yi;
+
+  int i = (int)yi - 1;
+  double *p[] = {
+    imgPixel(img, x, i+0),
+    imgPixel(img, x, i+1),
+    imgPixel(img, x, i+2),
+    imgPixel(img, x, i+3) };
+
+  pix[0] = cubicInterpolation(p[0][0], p[1][0], p[2][0], p[3][0], l);
+  pix[1] = cubicInterpolation(p[0][1], p[1][1], p[2][1], p[3][1], l);
+  pix[2] = cubicInterpolation(p[0][2], p[1][2], p[2][2], p[3][2], l);
+  pix[3] = cubicInterpolation(p[0][3], p[1][3], p[2][3], p[3][3], l);
+}
+
+
+double imgCubicCh(const Img *img, int ch, double x, double y) {
+  if (!img->data) return 0;
+  if (!(y > 0)) y = 0;
+  if (!(y < img->h)) y = img->h;
+  double yi = floor(y);
+  double l = y - yi;
+
+  int i = (int)yi - 1;
+  double p[] = {
+    imgCubicHCh(img, ch, x, i+0),
+    imgCubicHCh(img, ch, x, i+1),
+    imgCubicHCh(img, ch, x, i+2),
+    imgCubicHCh(img, ch, x, i+3) };
+
+  return cubicInterpolation(p[0], p[1], p[2], p[3], l);
+}
+
+
+void imgCubic(const Img *img, double x, double y, double *pix) {
+  if (!img->data) { pix[0] = pix[1] = pix[2] = pix[3] = 0; return; }
+  if (!(y > 0)) y = 0;
+  if (!(y < img->h)) y = img->h;
+  double yi = floor(y);
+  double l = y - yi;
+
+  int i = (int)yi - 1;
+  double p[4][4];
+  imgCubicH(img, x, i+0, p[0]);
+  imgCubicH(img, x, i+1, p[1]);
+  imgCubicH(img, x, i+2, p[2]);
+  imgCubicH(img, x, i+3, p[3]);
+
+  pix[0] = cubicInterpolation(p[0][0], p[1][0], p[2][0], p[3][0], l);
+  pix[1] = cubicInterpolation(p[0][1], p[1][1], p[2][1], p[3][1], l);
+  pix[2] = cubicInterpolation(p[0][2], p[1][2], p[2][2], p[3][2], l);
+  pix[3] = cubicInterpolation(p[0][3], p[1][3], p[2][3], p[3][3], l);
+}
+
+
+void imgFromInt(Img *img, int w, int h, const unsigned char *data) {
+  imgInit(img, w, h);
+  if (!img->data || !data) return;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; ++p)
+    *p = (*data++)/255.0;
+}
+
+
+unsigned char* imgToInt(Img *img) {
+  if (!img->data) return NULL;
+  unsigned char *data = (unsigned char*)malloc(sizeof(unsigned char)*4*img->w*img->h), *dp = data;
+  for(double *p = img->data, *e = p + 4*img->w*img->h; p < e; ++p)
+    *dp++ = *p > 0 ? (*p < 1 ? (unsigned char)floor(*p*255.9999999) : 255) : 0;
+  return data;
+}
+
diff --git a/img.ldr.inc.c b/img.ldr.inc.c
new file mode 100644
index 0000000..16025f3
--- /dev/null
+++ b/img.ldr.inc.c
@@ -0,0 +1,98 @@
+
+#include "img.ldr.tga.inc.c"
+
+#ifdef PNG_LIBPNG_VER_STRING
+#include "img.ldr.png.inc.c"
+#endif
+
+
+typedef int (*ImageLoadFunc)(const char *path, int *outWidth, int *outHeight, unsigned char **pixels);
+typedef int (*ImageSaveFunc)(const char *path, int width, int height, const unsigned char *pixels);
+
+
+int checkExt(const char *path, const char *ext) {
+  int lp = strlen(path);
+  int le = strlen(ext);
+  if (lp <= le) return FALSE;
+  if (path[lp-le-1] != '.') return FALSE;
+  for(const char *a = path+lp-le, *b = ext; *b; ++a, ++b)
+    if (tolower(*a) != tolower(*b)) return FALSE;
+  return TRUE;
+}
+
+
+int imageLoadAuto(const char *path, int *outWidth, int *outHeight, unsigned char **pixels) {
+  struct {
+    const char *ext;
+    ImageLoadFunc func;
+  } formats[] = {
+    { "tga", &imageLoadTga },
+    #if defined(HELI_HELIANTUS_H)
+    { "png", &imageLoad },
+    #elif defined(PNG_LIBPNG_VER_STRING)
+    { "png", &imageLoadPng },
+    #endif
+  };
+
+  int cnt = sizeof(formats)/sizeof(*formats);
+  for(int i = 0; i < cnt; ++i)
+    if (checkExt(path, formats[i].ext))
+      return formats[i].func(path, outWidth, outHeight, pixels);
+
+  fprintf(stderr, "Cannot load this image format (supported formats:");
+  for(int i = 0; i < cnt; ++i)
+    fprintf(stderr, " %s", formats[i].ext);
+  fprintf(stderr, "): %s\n", path);
+  fflush(stderr);
+  return FALSE;
+}
+
+
+int imageSaveAuto(const char *path, int width, int height, const unsigned char *pixels) {
+  struct {
+    const char *ext;
+    ImageSaveFunc func;
+  } formats[] = {
+    { "tga", &imageSaveTga },
+    #if defined(HELI_HELIANTUS_H)
+    { "png", (ImageSaveFunc)&imageSave },
+    #elif defined(PNG_LIBPNG_VER_STRING)
+    { "png", &imageSavePng },
+    #endif
+  };
+
+  int cnt = sizeof(formats)/sizeof(*formats);
+  for(int i = 0; i < cnt; ++i)
+    if (checkExt(path, formats[i].ext))
+      return formats[i].func(path, width, height, pixels);
+
+  fprintf(stderr, "Cannot save this image format (supported formats:");
+  for(int i = 0; i < cnt; ++i)
+    fprintf(stderr, " %s", formats[i].ext);
+  fprintf(stderr, "): %s\n", path);
+  fflush(stderr);
+  return FALSE;
+}
+
+
+int imgLoad(Img *img, const char *path) {
+  imgDestroy(img);
+  int w, h;
+  unsigned char *data;
+  if (!imageLoadAuto(path, &w, &h, &data))
+    { free(data); return FALSE; }
+  imgFromInt(img, w, h, data);
+  free(data);
+  return TRUE;
+}
+
+
+int imgSave(Img *img, const char *path) {
+  unsigned char *data = imgToInt(img);
+  if (!imageSaveAuto(path, img->w, img->h, data))
+    { free(data); return FALSE; }
+  free(data);
+  return TRUE;
+}
+
+
diff --git a/img.ldr.png.inc.c b/img.ldr.png.inc.c
new file mode 100644
index 0000000..a94d54e
--- /dev/null
+++ b/img.ldr.png.inc.c
@@ -0,0 +1,107 @@
+
+
+int imageLoadPng(const char *path, int *outWidth, int *outHeight, unsigned char **pixels) {
+  png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, 0, 0, 0);
+  png_infop info = png ? png_create_info_struct(png) : NULL;
+  if (!png) {
+    fprintf(stderr, "Cannot initialize PNG library\n");
+    fflush(stderr);
+    return FALSE;
+  }
+
+  FILE *f = fopen(path, "rb");
+  if (!f) {
+    fprintf(stderr, "Cannot open for read: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+  }
+
+  unsigned char header[8];
+  fread(header, sizeof(header), 1, f);
+  if (png_sig_cmp(header, 0, sizeof(header))) {
+    fclose(f);
+    fprintf(stderr, "Seems it is not PNG image: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+  }
+
+  png_init_io(png, f);
+  png_set_sig_bytes(png, sizeof(header));
+  png_read_info(png, info);
+
+  unsigned int width, height;
+  int bitdepth, colortype;
+  png_get_IHDR(png,info, &width, &height, &bitdepth, &colortype, 0, 0, 0);
+  if (colortype == PNG_COLOR_TYPE_PALETTE)
+    png_set_palette_to_rgb(png);
+  if (colortype == PNG_COLOR_TYPE_GRAY && bitdepth < 8)
+    png_set_expand_gray_1_2_4_to_8(png);
+  if (colortype == PNG_COLOR_TYPE_GRAY || colortype == PNG_COLOR_TYPE_GRAY_ALPHA)
+    png_set_gray_to_rgb(png);
+  if (png_get_valid(png, info, PNG_INFO_tRNS))
+    png_set_tRNS_to_alpha(png);
+  if (colortype == PNG_COLOR_TYPE_RGB)
+    png_set_filler(png, 0xff, PNG_FILLER_AFTER);
+  if (bitdepth == 16)
+    png_set_strip_16(png);
+  if (bitdepth < 8)
+    png_set_packing(png);
+
+  png_read_update_info(png, info);
+  png_get_IHDR(png, info, &width, &height, 0, 0, 0, 0, 0);
+
+  unsigned char *data = (unsigned char*)calloc(1, 4*width*height);
+  png_bytep *rows = (png_bytep*)malloc(height*sizeof(png_bytep));
+  for(unsigned int i = 0; i < height; i++)
+    rows[height-i-1] = data + i*4*width;
+  png_read_image(png, rows);
+  png_read_end(png, 0);
+  free(rows);
+
+  png_destroy_read_struct(&png, &info, 0);
+  fclose(f);
+
+  *outWidth = (int)width;
+  *outHeight = (int)height;
+  *pixels = data;
+  return TRUE;
+}
+
+
+int imageSavePng(const char *path, int width, int height, const unsigned char *pixels) {
+  png_structp png = png_create_write_struct(PNG_LIBPNG_VER_STRING, 0, 0, 0);
+  png_infop info = png ? png_create_info_struct(png) : NULL;
+  if (!png) {
+    fprintf(stderr, "Cannot initialize PNG library\n");
+    fflush(stderr);
+    return FALSE;
+  }
+
+  FILE *f = fopen(path, "wb");
+  if (!f) {
+    fprintf(stderr, "Cannot open for write: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+  }
+
+  png_init_io(png, f);
+  png_set_IHDR( png, info, width, height, 8, PNG_COLOR_TYPE_RGBA,
+                PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE );
+  png_write_info(png, info);
+  png_set_packing(png);
+
+  png_bytep *rows = (png_bytep*)malloc(height*sizeof(png_bytep));
+  for(int i = 0; i < height; i++)
+    rows[height-1-i] = (png_bytep)(pixels + i*4*width);
+  png_write_image(png, rows);
+  png_write_end(png, info);
+  free(rows);
+
+  png_destroy_write_struct(&png, &info);
+
+  fclose(f);
+  return TRUE;
+}
+
+
+
diff --git a/img.ldr.tga.inc.c b/img.ldr.tga.inc.c
new file mode 100644
index 0000000..0cd26b4
--- /dev/null
+++ b/img.ldr.tga.inc.c
@@ -0,0 +1,101 @@
+
+
+#pragma pack(push,1)
+typedef struct {
+  unsigned char  idLength;
+  unsigned char  colormapType;
+  unsigned char  imageType;
+  unsigned short colormapIndex;
+  unsigned short colormapLength;
+  unsigned char  colormapSize;
+  unsigned short xOrigin;
+  unsigned short yOrigin;
+  unsigned short width;
+  unsigned short height;
+  unsigned char  pixelSize;
+  unsigned char  attributes;
+} TgaHeader;
+#pragma pack(pop)
+
+
+int imageLoadTga(const char *path, int *outWidth, int *outHeight, unsigned char **pixels) {
+  *outWidth = *outHeight = 0; *pixels = NULL;
+
+	FILE *f = fopen(path, "rb");
+	if (!f) {
+    fprintf(stderr, "Cannot open for read: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+  }
+
+	TgaHeader header = {};
+	fread(&header, sizeof(header), 1, f); // todo: endianess
+
+	//printf("%d %d %d %d %d\n", header.imageType, header.width, header.height, header.colormapType, header.pixelSize);
+
+	if ( header.imageType != 2
+	  || !header.width
+	  || !header.height
+	  || header.colormapType
+	  || (header.pixelSize != 24 && header.pixelSize != 32)
+  ) {
+    fprintf(stderr, "Unsupported format. Only uncompressed 24 or 32 bits TGA are supported: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+	}
+
+  fseek(f, header.idLength, SEEK_CUR);
+
+  *outWidth = header.width;
+  *outHeight = header.height;
+  int rowSize = *outWidth*4;
+  int size = *outHeight*rowSize;
+  *pixels = (unsigned char*)calloc(1, size);
+  unsigned char *row = *pixels + size;
+  for(unsigned short r = header.height; r; --r, row -= rowSize) {
+    for(unsigned char *c = row - rowSize; c < row; c += 4) {
+      c[2] = fgetc(f);
+      c[1] = fgetc(f);
+      c[0] = fgetc(f);
+      c[3] = header.pixelSize == 32 ? fgetc(f) : 255;
+    }
+  }
+	fclose(f);
+
+	return TRUE;
+}
+
+
+int imageSaveTga(const char *path, int width, int height, const unsigned char *pixels) {
+  if (width <= 0 || height <= 0 || !pixels) return FALSE;
+
+	FILE *f = fopen(path, "wb");
+	if (!f) {
+    fprintf(stderr, "Cannot open file for write: %s\n", path);
+    fflush(stderr);
+    return FALSE;
+  }
+
+	TgaHeader header = {};
+	header.imageType = 2;
+	header.width = width;
+	header.height = height;
+	header.pixelSize = 32;
+	fwrite(&header, sizeof(header), 1, f); // todo: endianess
+
+  int rowSize = width*4;
+  int size = height*rowSize;
+  const unsigned char *row = pixels + size;
+  for(unsigned short r = header.height; r; --r, row -= rowSize) {
+    for(const unsigned char *c = row - rowSize; c < row; c += 4) {
+      fputc(c[2], f);
+      fputc(c[1], f);
+      fputc(c[0], f);
+      fputc(c[3], f);
+    }
+  }
+	fclose(f);
+
+	return TRUE;
+}
+