git.sesse.net Git - qscale/blob - qscale.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdlib.h>
   5
   6 #define CACHE_LINE_FACTOR 1
   7
   8 void myreadline(char *buf)
   9 {
  10         char *ptr;
  11
  12         if (fgets(buf, 256, stdin) == NULL) {
  13                 fprintf(stderr, "Premuture EOF\n");
  14                 exit(1);
  15         }
  16
  17         if (strlen(buf) == 0) {
  18                 return;
  19         }
  20
  21         ptr = &buf[strlen(buf) - 1];
  22         if (*ptr == '\n') {
  23                 *ptr-- = 0;
  24                 if (*ptr == '\r') {
  25                         *ptr = 0;
  26                 }
  27         }
  28 }
  29
  30 double sinc(double x)
  31 {
  32         // This is bad for very small x, should use power series instead.
  33         if (x == 0.0)
  34                 return 1.0;
  35         else
  36                 return sin(x) / x;
  37 }
  38
  39 double lanczos_tap(double x)
  40 {
  41         if (x < -3.0 || x > 3.0)
  42                 return 0.0;
  43         if (x < 0.0)
  44                 return sinc(-x*M_PI) * sinc(-x*M_PI / 3.0);
  45         else
  46                 return sinc(x*M_PI) * sinc(x*M_PI / 3.0);
  47 }
  48
  49
  50 struct pix_desc {
  51         unsigned start, end;
  52         unsigned startcoeff;
  53 };
  54
  55 void hscale(float *pix, unsigned char *npix, unsigned w, unsigned h, unsigned nw)
  56 {
  57         struct pix_desc *pd = (struct pix_desc *)malloc(nw * sizeof(struct pix_desc));
  58         int size_coeffs = 8;
  59         float *coeffs = (float *)malloc(size_coeffs * sizeof(float));
  60         int num_coeffs = 0;
  61         int x, y, sx;
  62         double sf = (double)w / (double)nw;
  63         double support = (w > nw) ? (3.0 * sf) : (3.0 / sf);
  64
  65         /* calculate the filter */
  66         for (x = 0; x < nw; ++x) {
  67                 int start = ceil(x * sf - support);
  68                 int end = floor(x * sf + support);
  69                 double sum = 0.0;
  70
  71                 if (start < 0) {
  72                         start = 0;
  73                 }
  74                 if (end > w - 1) {
  75                         end = w - 1;
  76                 }
  77
  78                 pd[x].start = start;
  79                 pd[x].end = end;
  80                 pd[x].startcoeff = num_coeffs;
  81
  82                 for (sx = start; sx <= end; ++sx) {
  83                         double nd = (w > nw) ? (sx/sf - x) : (sx - x*sf);
  84                         double f = lanczos_tap(nd);
  85                         if (num_coeffs == size_coeffs) {
  86                                 size_coeffs <<= 1;
  87                                 coeffs = (float *)realloc(coeffs, size_coeffs * sizeof(float));
  88                         }
  89
  90                         coeffs[num_coeffs++] = f;
  91                         sum += f;
  92                 }
  93
  94                 for (sx = start; sx <= end; ++sx) {
  95                         coeffs[pd[x].startcoeff + sx - start] /= sum;
  96                 }
  97         }
  98
  99         for (y = 0; y < h; ++y) {
 100                 float *sptr = pix + y*w;
 101                 unsigned char *dptr = npix + y*nw;
 102                 for (x = 0; x < nw; ++x) {
 103                         float acc = 0.0;
 104                         float *cf = &coeffs[pd[x].startcoeff];
 105                         unsigned sx;
 106                         unsigned char ch;
 107
 108                         for (sx = pd[x].start; sx <= pd[x].end; ++sx) {
 109                                 acc += sptr[sx] * *cf++;
 110                         }
 111
 112                         if (acc < 0.0)
 113                                 ch = 0;
 114                         else if (acc > 255.0)
 115                                 ch = 255;
 116                         else
 117                                 ch = (unsigned char)acc;
 118                         *dptr++ = ch;
 119                 }
 120         }
 121 }
 122
 123 void vscale(unsigned char *pix, float *npix, unsigned w, unsigned h, unsigned nh)
 124 {
 125         struct pix_desc *pd = (struct pix_desc *)malloc(nh * sizeof(struct pix_desc));
 126         int size_coeffs = 8;
 127         float *coeffs = (float *)malloc(size_coeffs * sizeof(float));
 128         int num_coeffs = 0;
 129         int x, y, sy;
 130         double sf = (double)h / (double)nh;
 131         double support = (h > nh) ? (3.0 * sf) : (3.0 / sf);
 132
 133         /* calculate the filter */
 134         for (y = 0; y < nh; ++y) {
 135                 int start = ceil(y * sf - support);
 136                 int end = floor(y * sf + support);
 137                 double sum = 0.0;
 138
 139                 if (start < 0) {
 140                         start = 0;
 141                 }
 142                 if (end > h - 1) {
 143                         end = h - 1;
 144                 }
 145
 146                 pd[y].start = start;
 147                 pd[y].end = end;
 148                 pd[y].startcoeff = num_coeffs;
 149
 150                 for (sy = start; sy <= end; ++sy) {
 151                         double nd = (h > nh) ? (sy/sf - y) : (sy - y*sf);
 152                         double f = lanczos_tap(nd);
 153                         if (num_coeffs == size_coeffs) {
 154                                 size_coeffs <<= 1;
 155                                 coeffs = (float *)realloc(coeffs, size_coeffs * sizeof(float));
 156                         }
 157
 158                         coeffs[num_coeffs++] = f;
 159                         sum += f;
 160                 }
 161
 162                 for (sy = start; sy <= end; ++sy) {
 163                         coeffs[pd[y].startcoeff + sy - start] /= sum;
 164                 }
 165         }
 166
 167 #if CACHE_LINE_FACTOR > 1
 168         for (x = 0; x < w; x += CACHE_LINE_FACTOR) {
 169                 unsigned char *sptr = pix + x;
 170                 float *dptr = npix + x;
 171                 for (y = 0; y < nh; ++y) {
 172                         int i;
 173                         float acc[CACHE_LINE_FACTOR];
 174                         for (i = 0; i < CACHE_LINE_FACTOR; ++i)
 175                                 acc[i] = 0.0;
 176                         float *cf = &coeffs[pd[y].startcoeff];
 177                         unsigned sy;
 178
 179                         for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
 180                                 for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
 181                                         acc[i] += sptr[sy * w + i] * *cf;
 182                                 }
 183                                 ++cf;
 184                         }
 185
 186                         for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
 187                                 dptr[i] = acc[i];
 188                         }
 189                         dptr += w;
 190                 }
 191         }
 192         for (x = (x/CACHE_LINE_FACTOR)*CACHE_LINE_FACTOR; x < w; ++x) {
 193 #else
 194         for (x = 0; x < w; ++x) {
 195 #endif
 196                 unsigned char *sptr = pix + x;
 197                 float *dptr = npix + x;
 198                 for (y = 0; y < nh; ++y) {
 199                         float acc = 0.0;
 200                         float *cf = &coeffs[pd[y].startcoeff];
 201                         unsigned sy;
 202
 203                         for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
 204                                 acc += sptr[sy * w] * *cf++;
 205                         }
 206
 207                         *dptr = acc;
 208                         dptr += w;
 209                 }
 210         }
 211 }
 212
 213 int main(void)
 214 {
 215         unsigned w, h;
 216         unsigned nw = 600, nh = 400;
 217         //unsigned nw = 3504, nh = 2336;
 218         char buf[256];
 219         unsigned char *pix;
 220         float *npix;
 221         unsigned char *npix2;
 222
 223         myreadline(buf);
 224         if (strcmp(buf, "P5") != 0) {
 225                 fprintf(stderr, "Error: Not a PGM file\n");
 226                 exit(1);
 227         }
 228         myreadline(buf);  /* comment */
 229         myreadline(buf);
 230
 231         if (sscanf(buf, "%u %u", &w, &h) != 2) {
 232                 fprintf(stderr, "Error: Expected width/height\n");
 233                 exit(1);
 234         }
 235
 236         myreadline(buf);
 237         if (strcmp(buf, "255") != 0) {
 238                 fprintf(stderr, "Error: Expected 8-bit\n");
 239                 exit(1);
 240         }
 241
 242         pix = (unsigned char *)malloc(w * h);
 243         if (pix == NULL) {
 244                 perror("malloc");
 245                 exit(1);
 246         }
 247
 248         if (fread(pix, w*h, 1, stdin) != 1) {
 249                 fprintf(stderr, "Error: Short read\n");
 250                 exit(1);
 251         }
 252
 253         npix = (float *)malloc(sizeof(float) * nw * h);
 254         if (npix == NULL) {
 255                 perror("malloc");
 256                 exit(1);
 257         }
 258         npix2 = (unsigned char *)malloc(nw * nh);
 259         if (npix2 == NULL) {
 260                 perror("malloc");
 261                 exit(1);
 262         }
 263
 264         vscale(pix, npix, w, h, nh);
 265         hscale(npix, npix2, w, nh, nw);
 266
 267         printf("P5\n# COMMENT: Made by qscale\n%u %u\n255\n", nw, nh);
 268         {
 269                 int y, x;
 270                 for (y = 0; y < nh; ++y) {
 271                         for (x = 0; x < nw; ++x) {
 272                                 putchar(npix2[y * nw + x]);
 273                         }
 274                 }
 275         }
 276
 277         return 0;
 278 }