blend.c 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <math.h>
  5. typedef struct {
  6. const char *fn;
  7. int xs;
  8. int ys;
  9. unsigned char *buf;
  10. } pgm;
  11. typedef struct {
  12. int xs;
  13. int ys;
  14. int *sum;
  15. int *count;
  16. } blendbuf;
  17. static void
  18. die (char *why)
  19. {
  20. fprintf (stderr, "%s\n", why);
  21. exit (1);
  22. }
  23. #define MAX_SIZE 65536
  24. pgm *load_pgm(const char *fn)
  25. {
  26. FILE *fi = fopen(fn, "rb");
  27. pgm *result;
  28. char buf[256];
  29. int xs, ys;
  30. int depth;
  31. if (fi == NULL)
  32. return NULL;
  33. fgets (buf, sizeof(buf), fi);
  34. if (buf[0] != 'P' || buf[1] != '5')
  35. die ("Need pgmraw image on input");
  36. xs = ys = 0;
  37. do
  38. fgets (buf, sizeof(buf), fi);
  39. while (buf[0] == '#');
  40. sscanf (buf, "%d %d", &xs, &ys);
  41. if (xs <= 0 || ys <= 0 || xs > MAX_SIZE || ys > MAX_SIZE)
  42. die ("Input image size out of range");
  43. do
  44. fgets (buf, sizeof(buf), fi);
  45. while (buf[0] == '#');
  46. sscanf (buf, "%d", &depth);
  47. if (depth != 255)
  48. die ("Only works with depth=255 images");
  49. result = (pgm *)malloc(sizeof(pgm));
  50. result->fn = fn;
  51. result->xs = xs;
  52. result->ys = ys;
  53. result->buf = (unsigned char *)malloc(xs * ys);
  54. fread(result->buf, 1, xs * ys, fi);
  55. fprintf(stderr, "loaded file %s %dx%d\n", fn, xs, ys);
  56. fclose(fi);
  57. return result;
  58. }
  59. int
  60. align_pgms(const pgm *p1, const pgm *p2, int *px, int *py)
  61. {
  62. int xo, yo;
  63. int xa, ya;
  64. int best = 0x7fffffff;
  65. xa = (p1->xs < p2->xs ? p1->xs : p2->xs) - 20;
  66. ya = (p1->ys < p2->ys ? p1->ys : p2->ys) - 20;
  67. for (yo = -10; yo <= 10; yo++)
  68. for (xo = -10; xo <= 10; xo++) {
  69. int sum = 0;
  70. int i, j;
  71. for (j = 0; j < ya; j++)
  72. for (i = 0; i < xa; i++) {
  73. int g1 = p1->buf[(j + 10) * p1->xs + i + 10];
  74. int g2 = p2->buf[(j + 10 - yo) * p2->xs + i - xo + 10];
  75. sum += (g1 - g2) * (g1 - g2);
  76. }
  77. if (sum < best) {
  78. best = sum;
  79. *px = xo;
  80. *py = yo;
  81. }
  82. }
  83. return best;
  84. }
  85. blendbuf *
  86. new_blendbuf(int xs, int ys)
  87. {
  88. blendbuf *result = (blendbuf *)malloc(sizeof(blendbuf));
  89. int i;
  90. result->xs = xs;
  91. result->ys = ys;
  92. result->sum = (int *)malloc(sizeof(int) * xs * ys);
  93. result->count = (int *)malloc(sizeof(int) * xs * ys);
  94. for (i = 0; i < xs * ys; i++) {
  95. result->sum[i] = 0;
  96. result->count[i] = 0;
  97. }
  98. return result;
  99. }
  100. void
  101. add_pgm(blendbuf *bb, pgm *p, int xo, int yo)
  102. {
  103. int i, j;
  104. for (j = 0; j < p->ys; j++) {
  105. if (j + yo >= 0 && j + yo < bb->ys) {
  106. for (i = 0; i < p->xs; i++) {
  107. if (i + xo >= 0 && i + xo < bb->xs) {
  108. int ix = (j + yo) * bb->xs + i + xo;
  109. bb->sum[ix] += p->buf[j * p->xs + i];
  110. bb->count[ix]++;
  111. }
  112. }
  113. }
  114. }
  115. }
  116. pgm *
  117. pgm_from_blendbuf(blendbuf *bb)
  118. {
  119. int xs = bb->xs;
  120. int ys = bb->ys;
  121. pgm *result = (pgm *)malloc(sizeof(pgm));
  122. int i, j;
  123. result->xs = xs;
  124. result->ys = ys;
  125. result->buf = (unsigned char *)malloc(xs * ys);
  126. for (j = 0; j < ys; j++) {
  127. for (i = 0; i < xs; i++) {
  128. int ix = j * xs + i;
  129. unsigned char g;
  130. if (bb->count[ix])
  131. g = (bb->sum[ix] + (bb->count[ix] >> 1)) / bb->count[ix];
  132. else
  133. g = 255;
  134. result->buf[ix] = g;
  135. }
  136. }
  137. return result;
  138. }
  139. pgm *
  140. pgm_from_blendbuf_enhanced(blendbuf *bb, double sharp, double gamma)
  141. {
  142. int xs = bb->xs;
  143. int ys = bb->ys;
  144. pgm *result = (pgm *)malloc(sizeof(pgm));
  145. int i, j;
  146. double ming = 255, maxg = 0;
  147. double *tmpbuf;
  148. result->xs = xs;
  149. result->ys = ys;
  150. result->buf = (unsigned char *)malloc(xs * ys);
  151. tmpbuf = (double *)malloc(xs * ys * sizeof(double));
  152. for (j = 0; j < ys; j++) {
  153. for (i = 0; i < xs; i++) {
  154. int ix = j * xs + i;
  155. double g;
  156. if (bb->count[ix]) {
  157. g = ((double)bb->sum[ix]) / bb->count[ix];
  158. if (g < ming) ming = g;
  159. if (g > maxg) maxg = g;
  160. } else
  161. g = 255.0;
  162. tmpbuf[ix] = g;
  163. }
  164. }
  165. for (j = 0; j < ys; j++) {
  166. for (i = 0; i < xs; i++) {
  167. int ix = j * xs + i;
  168. double g;
  169. if (bb->count[ix]) {
  170. int u, v;
  171. int cnt = 0;
  172. double sum = 0;
  173. for (v = -1; v <= 1; v++) {
  174. for (u = -1; u <= 1; u++) {
  175. if (i + u >= 0 && i + u < xs &&
  176. j + v >= 0 && j + v < ys &&
  177. bb->count[(j + v) * xs + i + u]) {
  178. sum += tmpbuf[(j + v) * xs + i + u];
  179. cnt += 1;
  180. }
  181. }
  182. }
  183. g = (1 + sharp) * tmpbuf[ix] - sharp * sum / cnt;
  184. g = (g - ming) / (maxg - ming);
  185. if (g < 0) g = 0;
  186. if (g > 1) g = 1;
  187. g = pow(g, gamma);
  188. } else
  189. g = 1;
  190. result->buf[ix] = (int)(255 * g + 0.5);
  191. }
  192. }
  193. free(tmpbuf);
  194. return result;
  195. }
  196. void
  197. print_pgm(pgm *p, FILE *f)
  198. {
  199. fprintf(f, "P5\n%d %d\n255\n", p->xs, p->ys);
  200. fwrite(p->buf, 1, p->xs * p->ys, f);
  201. }
  202. void
  203. threshold(pgm *p)
  204. {
  205. int i;
  206. for (i = 0; i < p->xs * p->ys; i++)
  207. p->buf[i] = p->buf[i] > 128 ? 1 : 0;
  208. }
  209. int
  210. classify(pgm **pgmlist, int n_pgm)
  211. {
  212. int *class = (int *)malloc(sizeof(int) * n_pgm);
  213. int n_class = 0;
  214. int i, j;
  215. int tshift = 4;
  216. for (i = 0; i < n_pgm; i++)
  217. class[i] = -1;
  218. for (i = 0; i < n_pgm; i++) {
  219. pgm *pi = pgmlist[i];
  220. if (class[i] == -1) {
  221. class[i] = n_class++;
  222. for (j = i + 1; j < n_pgm; j++) {
  223. pgm *pj = pgmlist[j];
  224. int xo, yo;
  225. int score;
  226. if (abs(pi->xs - pj->xs) < 10 &&
  227. abs(pi->ys - pj->ys) < 10) {
  228. score = align_pgms(pi, pj, &xo, &yo);
  229. if (score < ((pi->xs - 20) * (pi->ys - 20)) >> tshift) {
  230. class[j] = class[i];
  231. }
  232. }
  233. }
  234. }
  235. printf("%s: class%d\n", pi->fn, class[i]);
  236. fflush(stdout);
  237. }
  238. free(class);
  239. return 0;
  240. }
  241. static int
  242. intcompar(const void *a, const void *b) {
  243. return *((int *)a) - *((int *)b);
  244. }
  245. int
  246. do_blend(pgm **pgmlist, int n_pgm, int n_passes, float thresh)
  247. {
  248. blendbuf *bb;
  249. int pass;
  250. int i;
  251. pgm *base = pgmlist[0];
  252. int *scores, *scores2, *xos, *yos;
  253. scores = (int *)malloc(n_pgm * sizeof(int));
  254. scores2 = (int *)malloc(n_pgm * sizeof(int));
  255. xos = (int *)malloc(n_pgm * sizeof(int));
  256. yos = (int *)malloc(n_pgm * sizeof(int));
  257. for (pass = 0; pass < n_passes; pass++) {
  258. int scorethresh;
  259. bb = new_blendbuf(base->xs, base->ys);
  260. for (i = 0; i < n_pgm; i++) {
  261. int xo, yo;
  262. int score = align_pgms(base, pgmlist[i], &xo, &yo);
  263. fprintf(stderr, "%s: score = %d, offset = %d, %d\n",
  264. pgmlist[i]->fn, score, xo, yo);
  265. scores[i] = score;
  266. xos[i] = xo;
  267. yos[i] = yo;
  268. }
  269. if (pass == 0) {
  270. scorethresh = 0x7fffffff;
  271. } else {
  272. memcpy(scores2, scores, n_pgm * sizeof(int));
  273. qsort(scores2, n_pgm, sizeof(int), intcompar);
  274. scorethresh = scores2[(int)(thresh * n_pgm)];
  275. }
  276. for (i = 0; i < n_pgm; i++) {
  277. if (pass > 0)
  278. fprintf(stderr, "%s: score = %d %s\n",
  279. pgmlist[i]->fn, scores[i],
  280. scores[i] <= scorethresh ? "ok" : "-");
  281. if (scores[i] <= scorethresh)
  282. add_pgm(bb, pgmlist[i], xos[i], yos[i]);
  283. }
  284. if (pass == n_passes - 1)
  285. base = pgm_from_blendbuf_enhanced(bb, 5, 0.5);
  286. else
  287. base = pgm_from_blendbuf(bb);
  288. if (pass != n_passes - 1)
  289. fprintf(stderr, "\n");
  290. }
  291. free(scores);
  292. free(xos);
  293. free(yos);
  294. print_pgm(base, stdout);
  295. return 0;
  296. }
  297. int
  298. main(int argc, char **argv)
  299. {
  300. int i;
  301. int n_pgm = 0;
  302. pgm **pgmlist = (pgm **)malloc(sizeof(pgm *) * argc - 1);
  303. int do_class = 0;
  304. int n_pass = 2;
  305. float thresh = 0.90;
  306. for (i = 1; i < argc; i++) {
  307. if (!strcmp(argv[i], "-c"))
  308. do_class = 1;
  309. else
  310. pgmlist[n_pgm++] = load_pgm(argv[i]);
  311. }
  312. if (do_class) {
  313. for (i = 0; i < n_pgm; i++)
  314. threshold(pgmlist[i]);
  315. return classify(pgmlist, n_pgm);
  316. } else {
  317. return do_blend(pgmlist, n_pgm, n_pass, thresh);
  318. }
  319. return 0;
  320. }