sys/src/cmd/resample.c

   1 #include <u.h>
   2 #include <libc.h>
   3 #include <draw.h>
   4 #include <memdraw.h>
   5
   6 #define K2 7    /* from -.7 to +.7 inclusive, meaning .2 into each adjacent pixel */
   7 #define NK (2*K2+1)
   8 double K[NK];
   9
  10 double
  11 fac(int L)
  12 {
  13         int i, f;
  14
  15         f = 1;
  16         for(i=L; i>1; --i)
  17                 f *= i;
  18         return f;
  19 }
  20
  21 /*
  22  * i0(x) is the modified Bessel function, Σ (x/2)^2L / (L!)²
  23  * There are faster ways to calculate this, but we precompute
  24  * into a table so let's keep it simple.
  25  */
  26 double
  27 i0(double x)
  28 {
  29         double v;
  30         int L;
  31
  32         v = 1.0;
  33         for(L=1; L<10; L++)
  34                 v += pow(x/2., 2*L)/pow(fac(L), 2);
  35         return v;
  36 }
  37
  38 double
  39 kaiser(double x, double τ, double α)
  40 {
  41         if(fabs(x) > τ)
  42                 return 0.;
  43         return i0(α*sqrt(1-(x*x/(τ*τ))))/i0(α);
  44 }
  45
  46 void
  47 usage(void)
  48 {
  49         fprint(2, "usage: resample [-x xsize] [-y ysize] [imagefile]\n");
  50         fprint(2, "\twhere size is an integer or a percentage in the form 25%%\n");
  51         exits("usage");
  52 }
  53
  54 int
  55 getint(char *s, int *percent)
  56 {
  57         if(s == nil)
  58                 usage();
  59         *percent = (s[strlen(s)-1] == '%');
  60         if(*s == '+')
  61                 return atoi(s+1);
  62         if(*s == '-')
  63                 return -atoi(s+1);
  64         return atoi(s);
  65 }
  66
  67 void
  68 resamplex(uchar *in, int off, int d, int inx, uchar *out, int outx)
  69 {
  70         int i, x, k;
  71         double X, xx, v, rat;
  72
  73
  74         rat = (double)inx/(double)outx;
  75         for(x=0; x<outx; x++){
  76                 if(inx == outx){
  77                         /* don't resample if size unchanged */
  78                         out[off+x*d] = in[off+x*d];
  79                         continue;
  80                 }
  81                 v = 0.0;
  82                 X = x*rat;
  83                 for(k=-K2; k<=K2; k++){
  84                         xx = X + rat*k/10.;
  85                         i = xx;
  86                         if(i < 0)
  87                                 i = 0;
  88                         if(i >= inx)
  89                                 i = inx-1;
  90                         v += in[off+i*d] * K[K2+k];
  91                 }
  92                 out[off+x*d] = v;
  93         }
  94 }
  95
  96 void
  97 resampley(uchar **in, int off, int iny, uchar **out, int outy)
  98 {
  99         int y, i, k;
 100         double Y, yy, v, rat;
 101
 102         rat = (double)iny/(double)outy;
 103         for(y=0; y<outy; y++){
 104                 if(iny == outy){
 105                         /* don't resample if size unchanged */
 106                         out[y][off] = in[y][off];
 107                         continue;
 108                 }
 109                 v = 0.0;
 110                 Y = y*rat;
 111                 for(k=-K2; k<=K2; k++){
 112                         yy = Y + rat*k/10.;
 113                         i = yy;
 114                         if(i < 0)
 115                                 i = 0;
 116                         if(i >= iny)
 117                                 i = iny-1;
 118                         v += in[i][off] * K[K2+k];
 119                 }
 120                 out[y][off] = v;
 121         }
 122
 123 }
 124
 125 int
 126 max(int a, int b)
 127 {
 128         if(a > b)
 129                 return a;
 130         return b;
 131 }
 132
 133 Memimage*
 134 resample(int xsize, int ysize, Memimage *m)
 135 {
 136         int i, j, bpl, nchan;
 137         Memimage *new;
 138         uchar **oscan, **nscan;
 139
 140         new = allocmemimage(Rect(0, 0, xsize, ysize), m->chan);
 141         if(new == nil)
 142                 sysfatal("can't allocate new image: %r");
 143
 144         oscan = malloc(Dy(m->r)*sizeof(uchar*));
 145         nscan = malloc(max(ysize, Dy(m->r))*sizeof(uchar*));
 146         if(oscan == nil || nscan == nil)
 147                 sysfatal("can't allocate: %r");
 148
 149         /* unload original image into scan lines */
 150         bpl = bytesperline(m->r, m->depth);
 151         for(i=0; i<Dy(m->r); i++){
 152                 oscan[i] = malloc(bpl);
 153                 if(oscan[i] == nil)
 154                         sysfatal("can't allocate: %r");
 155                 j = unloadmemimage(m, Rect(m->r.min.x, m->r.min.y+i, m->r.max.x, m->r.min.y+i+1), oscan[i], bpl);
 156                 if(j != bpl)
 157                         sysfatal("unloadmemimage");
 158         }
 159
 160         /* allocate scan lines for destination. we do y first, so need at least Dy(m->r) lines */
 161         bpl = bytesperline(Rect(0, 0, xsize, Dy(m->r)), m->depth);
 162         for(i=0; i<max(ysize, Dy(m->r)); i++){
 163                 nscan[i] = malloc(bpl);
 164                 if(nscan[i] == nil)
 165                         sysfatal("can't allocate: %r");
 166         }
 167
 168         /* resample in X */
 169         nchan = m->depth/8;
 170         for(i=0; i<Dy(m->r); i++){
 171                 for(j=0; j<nchan; j++){
 172                         if(j==0 && m->chan==XRGB32)
 173                                 continue;
 174                         resamplex(oscan[i], j, nchan, Dx(m->r), nscan[i], xsize);
 175                 }
 176                 free(oscan[i]);
 177                 oscan[i] = nscan[i];
 178                 nscan[i] = malloc(bpl);
 179                 if(nscan[i] == nil)
 180                         sysfatal("can't allocate: %r");
 181         }
 182
 183         /* resample in Y */
 184         for(i=0; i<xsize; i++)
 185                 for(j=0; j<nchan; j++)
 186                         resampley(oscan, nchan*i+j, Dy(m->r), nscan, ysize);
 187
 188         /* pack data into destination */
 189         bpl = bytesperline(new->r, m->depth);
 190         for(i=0; i<ysize; i++){
 191                 j = loadmemimage(new, Rect(0, i, xsize, i+1), nscan[i], bpl);
 192                 if(j != bpl)
 193                         sysfatal("loadmemimage: %r");
 194         }
 195         return new;
 196 }
 197
 198 void
 199 main(int argc, char *argv[])
 200 {
 201         int i, fd, xsize, ysize, xpercent, ypercent;
 202         Rectangle rparam;
 203         Memimage *m, *new, *t1, *t2;
 204         char *file;
 205         ulong tchan;
 206         char tmp[100];
 207         double v;
 208
 209         for(i=-K2; i<=K2; i++){
 210                 K[K2+i] = kaiser(i/10., K2/10., 4.);
 211 //              print("%g %g\n", i/10., K[K2+i]);
 212         }
 213
 214         /* normalize */
 215         v = 0.0;
 216         for(i=0; i<NK; i++)
 217                 v += K[i];
 218         for(i=0; i<NK; i++)
 219                 K[i] /= v;
 220
 221         memimageinit();
 222         memset(&rparam, 0, sizeof rparam);
 223         xsize = ysize = 0;
 224         xpercent = ypercent = 0;
 225
 226         ARGBEGIN{
 227         case 'a':       /* compatibility; equivalent to just -x or -y */
 228                 if(xsize != 0 || ysize != 0)
 229                         usage();
 230                 xsize = getint(ARGF(), &xpercent);
 231                 if(xsize <= 0)
 232                         usage();
 233                 ysize = xsize;
 234                 ypercent = xpercent;
 235                 break;
 236         case 'x':
 237                 if(xsize != 0)
 238                         usage();
 239                 xsize = getint(ARGF(), &xpercent);
 240                 if(xsize <= 0)
 241                         usage();
 242                 break;
 243         case 'y':
 244                 if(ysize != 0)
 245                         usage();
 246                 ysize = getint(ARGF(), &ypercent);
 247                 if(ysize <= 0)
 248                         usage();
 249                 break;
 250         default:
 251                 usage();
 252         }ARGEND
 253
 254         if(xsize == 0 && ysize == 0)
 255                 usage();
 256
 257         file = "<stdin>";
 258         fd = 0;
 259         if(argc > 1)
 260                 usage();
 261         else if(argc == 1){
 262                 file = argv[0];
 263                 fd = open(file, OREAD);
 264                 if(fd < 0)
 265                         sysfatal("can't open %s: %r", file);
 266         }
 267
 268         m = readmemimage(fd);
 269         if(m == nil)
 270                 sysfatal("can't read %s: %r", file);
 271
 272         if(xpercent)
 273                 xsize = Dx(m->r)*xsize/100;
 274         if(ypercent)
 275                 ysize = Dy(m->r)*ysize/100;
 276         if(ysize == 0)
 277                 ysize = (xsize * Dy(m->r)) / Dx(m->r);
 278         if(xsize == 0)
 279                 xsize = (ysize * Dx(m->r)) / Dy(m->r);
 280
 281         new = nil;
 282         switch(m->chan){
 283
 284         case GREY8:
 285         case RGB24:
 286         case RGBA32:
 287         case ARGB32:
 288         case XRGB32:
 289                 new = resample(xsize, ysize, m);
 290                 break;
 291
 292         case CMAP8:
 293         case RGB15:
 294         case RGB16:
 295                 tchan = RGB24;
 296                 goto Convert;
 297
 298         case GREY1:
 299         case GREY2:
 300         case GREY4:
 301                 tchan = GREY8;
 302         Convert:
 303                 /* use library to convert to byte-per-chan form, then convert back */
 304                 t1 = allocmemimage(m->r, tchan);
 305                 if(t1 == nil)
 306                         sysfatal("can't allocate temporary image: %r");
 307                 memimagedraw(t1, t1->r, m, m->r.min, nil, ZP, S);
 308                 t2 = resample(xsize, ysize, t1);
 309                 freememimage(t1);
 310                 new = allocmemimage(Rect(0, 0, xsize, ysize), m->chan);
 311                 if(new == nil)
 312                         sysfatal("can't allocate new image: %r");
 313                 /* should do error diffusion here */
 314                 memimagedraw(new, new->r, t2, t2->r.min, nil, ZP, S);
 315                 freememimage(t2);
 316                 break;
 317
 318         default:
 319                 sysfatal("can't handle channel type %s", chantostr(tmp, m->chan));
 320         }
 321
 322         assert(new);
 323         if(writememimage(1, new) < 0)
 324                 sysfatal("write error on output: %r");
 325
 326         exits(nil);
 327 }