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 }