tresample.c - plan9port - [fork] Plan 9 from user space
HTML git clone git://src.adamsgaard.dk/plan9port
DIR Log
DIR Files
DIR Refs
DIR README
DIR LICENSE
---
tresample.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 }