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