--- /dev/null
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+#define CACHE_LINE_FACTOR 1
+
+void myreadline(char *buf)
+{
+ char *ptr;
+
+ if (fgets(buf, 256, stdin) == NULL) {
+ fprintf(stderr, "Premuture EOF\n");
+ exit(1);
+ }
+
+ if (strlen(buf) == 0) {
+ return;
+ }
+
+ ptr = &buf[strlen(buf) - 1];
+ if (*ptr == '\n') {
+ *ptr-- = 0;
+ if (*ptr == '\r') {
+ *ptr = 0;
+ }
+ }
+}
+
+double sinc(double x)
+{
+ // This is bad for very small x, should use power series instead.
+ if (x == 0.0)
+ return 1.0;
+ else
+ return sin(x) / x;
+}
+
+double lanczos_tap(double x)
+{
+ if (x < -3.0 || x > 3.0)
+ return 0.0;
+ if (x < 0.0)
+ return sinc(-x*M_PI) * sinc(-x*M_PI / 3.0);
+ else
+ return sinc(x*M_PI) * sinc(x*M_PI / 3.0);
+}
+
+
+struct pix_desc {
+ unsigned start, end;
+ unsigned startcoeff;
+};
+
+void hscale(unsigned char *pix, float *npix, unsigned w, unsigned h, unsigned nw)
+{
+ struct pix_desc *pd = (struct pix_desc *)malloc(nw * sizeof(struct pix_desc));
+ int size_coeffs = 8;
+ float *coeffs = (float *)malloc(size_coeffs * sizeof(float));
+ int num_coeffs = 0;
+ int x, y, sx;
+ double sf = (double)w / (double)nw;
+ double support = 3.0 * sf;
+
+ /* calculate the filter */
+ for (x = 0; x < nw; ++x) {
+ int start = ceil(x * sf - support);
+ int end = floor(x * sf + support);
+ double sum = 0.0;
+
+ if (start < 0) {
+ start = 0;
+ }
+ if (end > w - 1) {
+ end = w - 1;
+ }
+
+ pd[x].start = start;
+ pd[x].end = end;
+ pd[x].startcoeff = num_coeffs;
+
+ for (sx = start; sx <= end; ++sx) {
+ double nd = sx/sf - x;
+ double f = lanczos_tap(nd);
+ if (num_coeffs == size_coeffs) {
+ size_coeffs <<= 1;
+ coeffs = (float *)realloc(coeffs, size_coeffs * sizeof(float));
+ }
+
+ coeffs[num_coeffs++] = f;
+ sum += f;
+ }
+
+ for (sx = start; sx <= end; ++sx) {
+ coeffs[pd[x].startcoeff + sx - start] /= sum;
+ }
+ }
+
+ for (y = 0; y < h; ++y) {
+ unsigned char *sptr = pix + y*w;
+ float *dptr = npix + y*nw;
+ for (x = 0; x < nw; ++x) {
+ float acc = 0.0;
+ float *cf = &coeffs[pd[x].startcoeff];
+ unsigned sx;
+
+ for (sx = pd[x].start; sx <= pd[x].end; ++sx) {
+ acc += sptr[sx] * *cf++;
+ }
+ *dptr++ = acc;
+ }
+ }
+}
+
+void vscale(float *pix, unsigned char *npix, unsigned w, unsigned h, unsigned nh)
+{
+ struct pix_desc *pd = (struct pix_desc *)malloc(nh * sizeof(struct pix_desc));
+ int size_coeffs = 8;
+ float *coeffs = (float *)malloc(size_coeffs * sizeof(float));
+ int num_coeffs = 0;
+ int x, y, sy;
+ double sf = (double)h / (double)nh;
+ double support = 3.0 * sf;
+
+ /* calculate the filter */
+ for (y = 0; y < nh; ++y) {
+ int start = ceil(y * sf - support);
+ int end = floor(y * sf + support);
+ double sum = 0.0;
+
+ if (start < 0) {
+ start = 0;
+ }
+ if (end > h - 1) {
+ end = h - 1;
+ }
+
+ pd[y].start = start;
+ pd[y].end = end;
+ pd[y].startcoeff = num_coeffs;
+
+ for (sy = start; sy <= end; ++sy) {
+ double nd = sy/sf - y;
+ double f = lanczos_tap(nd);
+ if (num_coeffs == size_coeffs) {
+ size_coeffs <<= 1;
+ coeffs = (float *)realloc(coeffs, size_coeffs * sizeof(float));
+ }
+
+ coeffs[num_coeffs++] = f;
+ sum += f;
+ }
+
+ for (sy = start; sy <= end; ++sy) {
+ coeffs[pd[y].startcoeff + sy - start] /= sum;
+ }
+ }
+#if CACHE_LINE_FACTOR > 1
+ for (x = 0; x < w; x += CACHE_LINE_FACTOR) {
+ float *sptr = pix + x;
+ unsigned char *dptr = npix + x;
+ for (y = 0; y < nh; ++y) {
+ int i;
+ float acc[CACHE_LINE_FACTOR];
+ for (i = 0; i < CACHE_LINE_FACTOR; ++i)
+ acc[i] = 0.0;
+ float *cf = &coeffs[pd[y].startcoeff];
+ unsigned sy;
+ unsigned char ch;
+
+ for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
+ for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
+ acc[i] += sptr[sy * w + i] * *cf;
+ }
+ ++cf;
+ }
+
+ for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
+ if (acc[i] < 0.0)
+ ch = 0;
+ else if (acc[i] > 255.0)
+ ch = 255;
+ else
+ ch = (unsigned char)acc[i];
+ dptr[i] = ch;
+ }
+ dptr += w;
+ }
+ }
+#endif
+ for (x = (x/CACHE_LINE_FACTOR)*CACHE_LINE_FACTOR; x < w; ++x) {
+ float *sptr = pix + x;
+ unsigned char *dptr = npix + x;
+ for (y = 0; y < nh; ++y) {
+ float acc = 0.0;
+ float *cf = &coeffs[pd[y].startcoeff];
+ unsigned sy;
+ unsigned char ch;
+
+ for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
+ acc += sptr[sy * w] * *cf++;
+ }
+
+ if (acc < 0.0)
+ ch = 0;
+ else if (acc > 255.0)
+ ch = 255;
+ else
+ ch = (unsigned char)acc;
+
+ *dptr = ch;
+ dptr += w;
+ }
+ }
+}
+
+int main(void)
+{
+ unsigned w, h;
+ unsigned nw = 600, nh = 400;
+ char buf[256];
+ unsigned char *pix;
+ float *npix;
+ unsigned char *npix2;
+
+ myreadline(buf);
+ if (strcmp(buf, "P5") != 0) {
+ fprintf(stderr, "Error: Not a PGM file\n");
+ exit(1);
+ }
+ myreadline(buf); /* comment */
+ myreadline(buf);
+
+ if (sscanf(buf, "%u %u", &w, &h) != 2) {
+ fprintf(stderr, "Error: Expected width/height\n");
+ exit(1);
+ }
+
+ myreadline(buf);
+ if (strcmp(buf, "255") != 0) {
+ fprintf(stderr, "Error: Expected 8-bit\n");
+ exit(1);
+ }
+
+ pix = (unsigned char *)malloc(w * h);
+ if (pix == NULL) {
+ perror("malloc");
+ exit(1);
+ }
+
+ if (fread(pix, w*h, 1, stdin) != 1) {
+ fprintf(stderr, "Error: Short read\n");
+ exit(1);
+ }
+
+ npix = (float *)malloc(sizeof(float) * nw * h);
+ if (npix == NULL) {
+ perror("malloc");
+ exit(1);
+ }
+ npix2 = (unsigned char *)malloc(nw * nh);
+ if (npix2 == NULL) {
+ perror("malloc");
+ exit(1);
+ }
+
+ hscale(pix, npix, w, h, nw);
+ vscale(npix, npix2, nw, h, nh);
+
+ printf("P5\n# COMMENT: Made by qscale\n%u %u\n255\n", nw, nh);
+ {
+ int y, x;
+ for (y = 0; y < nh; ++y) {
+ for (x = 0; x < nw; ++x) {
+ putchar(npix2[y * nw + x]);
+ }
+ }
+ }
+
+ return 0;
+}
projects / qscale / commitdiff