resample.c (6283B)
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 tau, double alpha) 40 { 41 if(fabs(x) > tau) 42 return 0.; 43 return i0(alpha*sqrt(1-(x*x/(tau*tau))))/i0(alpha); 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 }