From: sgunderson@bigfoot.com <> Date: Sun, 3 Feb 2008 16:06:46 +0000 (+0100) Subject: Initial version. X-Git-Url: https://git.sesse.net/?p=qscale;a=commitdiff_plain;h=44d2d19e679c578800f420739385a159b8246949 Initial version. --- 44d2d19e679c578800f420739385a159b8246949 diff --git a/qscale.c b/qscale.c new file mode 100644 index 0000000..61f49c4 --- /dev/null +++ b/qscale.c @@ -0,0 +1,281 @@ +#include +#include +#include +#include + +#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; +}