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(unsigned char *pix, float *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 = 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 = sx/sf - x;
  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                 unsigned char *sptr = pix + y*w;
 101                 float *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
 107                         for (sx = pd[x].start; sx <= pd[x].end; ++sx) {
 108                                 acc += sptr[sx] * *cf++;
 109                         }
 110                         *dptr++ = acc;
 111                 }
 112         }
 113 }
 114
 115 void vscale(float *pix, unsigned char *npix, unsigned w, unsigned h, unsigned nh)
 116 {
 117         struct pix_desc *pd = (struct pix_desc *)malloc(nh * sizeof(struct pix_desc));
 118         int size_coeffs = 8;
 119         float *coeffs = (float *)malloc(size_coeffs * sizeof(float));
 120         int num_coeffs = 0;
 121         int x, y, sy;
 122         double sf = (double)h / (double)nh;
 123         double support = 3.0 * sf;
 124
 125         /* calculate the filter */
 126         for (y = 0; y < nh; ++y) {
 127                 int start = ceil(y * sf - support);
 128                 int end = floor(y * sf + support);
 129                 double sum = 0.0;
 130
 131                 if (start < 0) {
 132                         start = 0;
 133                 }
 134                 if (end > h - 1) {
 135                         end = h - 1;
 136                 }
 137
 138                 pd[y].start = start;
 139                 pd[y].end = end;
 140                 pd[y].startcoeff = num_coeffs;
 141
 142                 for (sy = start; sy <= end; ++sy) {
 143                         double nd = sy/sf - y;
 144                         double f = lanczos_tap(nd);
 145                         if (num_coeffs == size_coeffs) {
 146                                 size_coeffs <<= 1;
 147                                 coeffs = (float *)realloc(coeffs, size_coeffs * sizeof(float));
 148                         }
 149
 150                         coeffs[num_coeffs++] = f;
 151                         sum += f;
 152                 }
 153
 154                 for (sy = start; sy <= end; ++sy) {
 155                         coeffs[pd[y].startcoeff + sy - start] /= sum;
 156                 }
 157         }
 158 #if CACHE_LINE_FACTOR > 1
 159         for (x = 0; x < w; x += CACHE_LINE_FACTOR) {
 160                 float *sptr = pix + x;
 161                 unsigned char *dptr = npix + x;
 162                 for (y = 0; y < nh; ++y) {
 163                         int i;
 164                         float acc[CACHE_LINE_FACTOR];
 165                         for (i = 0; i < CACHE_LINE_FACTOR; ++i)
 166                                 acc[i] = 0.0;
 167                         float *cf = &coeffs[pd[y].startcoeff];
 168                         unsigned sy;
 169                         unsigned char ch;
 170
 171                         for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
 172                                 for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
 173                                         acc[i] += sptr[sy * w + i] * *cf;
 174                                 }
 175                                 ++cf;
 176                         }
 177
 178                         for (i = 0; i < CACHE_LINE_FACTOR; ++i) {
 179                                 if (acc[i] < 0.0)
 180                                         ch = 0;
 181                                 else if (acc[i] > 255.0)
 182                                         ch = 255;
 183                                 else
 184                                         ch = (unsigned char)acc[i];
 185                                 dptr[i] = ch;
 186                         }
 187                         dptr += w;
 188                 }
 189         }
 190 #endif
 191         for (x = (x/CACHE_LINE_FACTOR)*CACHE_LINE_FACTOR; x < w; ++x) {
 192                 float *sptr = pix + x;
 193                 unsigned char *dptr = npix + x;
 194                 for (y = 0; y < nh; ++y) {
 195                         float acc = 0.0;
 196                         float *cf = &coeffs[pd[y].startcoeff];
 197                         unsigned sy;
 198                         unsigned char ch;
 199
 200                         for (sy = pd[y].start; sy <= pd[y].end; ++sy) {
 201                                 acc += sptr[sy * w] * *cf++;
 202                         }
 203
 204                         if (acc < 0.0)
 205                                 ch = 0;
 206                         else if (acc > 255.0)
 207                                 ch = 255;
 208                         else
 209                                 ch = (unsigned char)acc;
 210
 211                         *dptr = ch;
 212                         dptr += w;
 213                 }
 214         }
 215 }
 216
 217 int main(void)
 218 {
 219         unsigned w, h;
 220         unsigned nw = 600, nh = 400;
 221         char buf[256];
 222         unsigned char *pix;
 223         float *npix;
 224         unsigned char *npix2;
 225
 226         myreadline(buf);
 227         if (strcmp(buf, "P5") != 0) {
 228                 fprintf(stderr, "Error: Not a PGM file\n");
 229                 exit(1);
 230         }
 231         myreadline(buf);  /* comment */
 232         myreadline(buf);
 233
 234         if (sscanf(buf, "%u %u", &w, &h) != 2) {
 235                 fprintf(stderr, "Error: Expected width/height\n");
 236                 exit(1);
 237         }
 238
 239         myreadline(buf);
 240         if (strcmp(buf, "255") != 0) {
 241                 fprintf(stderr, "Error: Expected 8-bit\n");
 242                 exit(1);
 243         }
 244
 245         pix = (unsigned char *)malloc(w * h);
 246         if (pix == NULL) {
 247                 perror("malloc");
 248                 exit(1);
 249         }
 250
 251         if (fread(pix, w*h, 1, stdin) != 1) {
 252                 fprintf(stderr, "Error: Short read\n");
 253                 exit(1);
 254         }
 255
 256         npix = (float *)malloc(sizeof(float) * nw * h);
 257         if (npix == NULL) {
 258                 perror("malloc");
 259                 exit(1);
 260         }
 261         npix2 = (unsigned char *)malloc(nw * nh);
 262         if (npix2 == NULL) {
 263                 perror("malloc");
 264                 exit(1);
 265         }
 266
 267         hscale(pix, npix, w, h, nw);
 268         vscale(npix, npix2, nw, h, nh);
 269
 270         printf("P5\n# COMMENT: Made by qscale\n%u %u\n255\n", nw, nh);
 271         {
 272                 int y, x;
 273                 for (y = 0; y < nh; ++y) {
 274                         for (x = 0; x < nw; ++x) {
 275                                 putchar(npix2[y * nw + x]);
 276                         }
 277                 }
 278         }
 279
 280         return 0;
 281 }