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