resample.c 6.1 KB

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