Mercurial > hg > wm
annotate Meerwald-dir/wm_bruyn_d.c @ 25:5a57a145bccb default tip
Added tag v0.7 for changeset 9f20bce6184e
author | Peter Meerwald-Stadler <pmeerw@pmeerw.net> |
---|---|
date | Fri, 20 Dec 2024 13:32:15 +0100 |
parents | 9f20bce6184e |
children |
rev | line source |
---|---|
0 | 1 #include "wm.h" |
2 #include "signature.h" | |
3 #include "coord.h" | |
4 #include "gray.h" | |
5 #include "sort.h" | |
6 #include "bruyn_common.h" | |
24
9f20bce6184e
move directories, support netpbm 11
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
20
diff
changeset
|
7 #include "netpbm/pgm.h" |
0 | 8 |
9 char *progname; | |
10 | |
11 // prints out program's parameters | |
12 void usage(void) { | |
13 fprintf(stderr, "usage: %s [-b n] [-h] [-k] [-n n] [-o file] [-pP n] [-q n] [-tT n] [-v n] -s file file\n", progname); | |
14 fprintf(stderr, "\t-b n\t\tblock size\n"); | |
15 fprintf(stderr, "\t-h\t\tprint usage\n"); | |
16 fprintf(stderr, "\t-k\t\tdisable block skipping\n"); | |
17 fprintf(stderr, "\t-n n\t\tnumber of signature bits to detect\n"); | |
18 fprintf(stderr, "\t-o file\t\textracted signature file\n"); | |
19 fprintf(stderr, "\t-p n\t\tpattern type for zone 1\n"); | |
20 fprintf(stderr, "\t-P n\t\tpattern type for zone 2\n"); | |
21 fprintf(stderr, "\t-q n\t\tsignature strength\n"); | |
22 fprintf(stderr, "\t-s file\t\tembedded signature\n"); | |
23 fprintf(stderr, "\t-t n\t\tthreshold for noise\n"); | |
24 fprintf(stderr, "\t-T n\t\tthreshold for slope\n"); | |
25 fprintf(stderr, "\t-v n\t\tverbosity level\n"); | |
26 exit(0); | |
27 } | |
28 | |
29 int main(int argc, char *argv[]) { | |
30 FILE *in = stdin; | |
31 FILE *out = stdout; | |
32 FILE *sig = NULL; | |
33 | |
34 gray** image; | |
35 gray **block; | |
36 gray **zone; | |
37 gray **category1, **category2; | |
38 gray maxval; | |
39 double *slope; | |
8 | 40 int rows, cols, format; |
0 | 41 int c; |
42 int i, j; | |
43 int n; | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
16
diff
changeset
|
44 int row; |
0 | 45 int n_block; |
46 | |
47 char signature_name[MAXPATHLEN]; | |
48 char input_name[MAXPATHLEN] = "(stdin)"; | |
49 char output_name[MAXPATHLEN] = "(stdout)"; | |
50 | |
51 double quality = 0.0; | |
52 double threshold_noise = 0.0; | |
53 double threshold_slope = 0.0; | |
54 int pattern1 = 0; | |
55 int pattern2 = 0; | |
56 int blocksize = 0; | |
57 int seed; | |
58 | |
59 int verbose = 0; | |
60 int skipping = 0; | |
61 | |
62 struct coords *coords; | |
63 | |
64 progname = argv[0]; | |
65 | |
66 pgm_init(&argc, argv); wm_init2(); | |
67 | |
68 // parse command line and set options | |
69 while ((c = getopt(argc, argv, "b:h?n:o:p:P:q:s:t:T:v:k")) != EOF) { | |
70 switch (c) { | |
71 case 'h': | |
72 case '?': | |
73 usage(); | |
74 break; | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
16
diff
changeset
|
75 case 'k': |
0 | 76 skipping = 1; |
77 break; | |
78 case 'n': | |
79 nbit_signature = atoi(optarg); | |
80 if (nbit_signature <= 0 || nbit_signature > NBITSIGNATURE) { | |
81 fprintf(stderr, "%s: invalid signature length %d\n", progname, nbit_signature); | |
82 exit(1); | |
83 } | |
84 break; | |
85 case 'o': | |
86 if ((out = fopen(optarg, "wb")) == NULL) { | |
87 fprintf(stderr, "%s: unable to open output file %s\n", progname, optarg); | |
88 exit(1); | |
89 } | |
90 strcpy(output_name, optarg); | |
91 break; | |
92 case 'p': | |
93 pattern1 = atoi(optarg); | |
94 if (pattern1 <= 0 || pattern1 > NPATTERN) { | |
95 fprintf(stderr, "%s: pattern type out of range\n", progname); | |
96 exit(1); | |
97 } | |
98 break; | |
99 case 'P': | |
100 pattern2 = atoi(optarg); | |
101 if (pattern2 <= 0 || pattern2 > 3) { | |
102 fprintf(stderr, "%s: pattern type out of range\n", progname); | |
103 exit(1); | |
104 } | |
105 break; | |
106 case 'q': | |
107 quality = atof(optarg); | |
108 if (quality <= 0) { | |
109 fprintf(stderr, "%s: quality factor %f out of range\n", progname, quality); | |
110 } | |
111 break; | |
112 case 's': | |
113 if ((sig = fopen(optarg, "r")) == NULL) { | |
114 fprintf(stderr, "%s: unable to open signature file %s\n", progname, optarg); | |
115 exit(1); | |
116 } | |
117 strcpy(signature_name, optarg); | |
118 break; | |
119 case 't': | |
120 threshold_noise = atof(optarg); | |
121 if (threshold_noise <= 0) { | |
122 fprintf(stderr, "%s: noise threshold %f out of range\n", progname, threshold_noise); | |
123 } | |
124 break; | |
125 case 'T': | |
126 threshold_slope = atof(optarg); | |
127 if (threshold_slope <= 0) { | |
128 fprintf(stderr, "%s: slope threshold %f out of range\n", progname, threshold_slope); | |
129 } | |
130 break; | |
131 case 'v': | |
132 verbose = atoi(optarg); | |
133 if (verbose < 0) { | |
134 fprintf(stderr, "%s: verbosity level %d out of range\n",progname, verbose); | |
135 exit(1); | |
136 } | |
137 break; | |
138 } | |
139 } | |
140 | |
141 argc -= optind; | |
142 argv += optind; | |
143 | |
144 if (argc > 1) { | |
145 usage(); | |
146 exit(1); | |
147 } | |
148 | |
149 // open input image file or read from stdin | |
8 | 150 if (argc == 1 && *argv[0] != '-') { |
0 | 151 if ((in = fopen(argv[0], "rb")) == NULL) { |
152 fprintf(stderr, "%s: unable to open input file %s\n", progname, argv[0]); | |
153 exit(1); | |
154 } | |
155 else | |
156 strcpy(input_name, argv[0]); | |
8 | 157 } |
0 | 158 |
159 // read signature file and set options | |
160 // command line options override signature file options | |
161 if (sig) { | |
8 | 162 char line[1024]; |
0 | 163 fgets(line, sizeof(line), sig); |
164 if (strspn(line, "BRSG") >= 4) { | |
165 if (nbit_signature == 0) | |
166 fscanf(sig, "%d\n", &nbit_signature); | |
167 else | |
168 fscanf(sig, "%*d\n"); | |
169 if (skipping == 0) | |
170 fscanf(sig, "%d\n", &skipping); | |
171 else | |
172 fscanf(sig, "%*d\n"); | |
173 if (pattern1 == 0) | |
174 fscanf(sig, "%d\n", &pattern1); | |
175 else | |
176 fscanf(sig, "%*d\n"); | |
177 if (pattern2 == 0) | |
178 fscanf(sig, "%d\n", &pattern2); | |
179 else | |
180 fscanf(sig, "%*d\n"); | |
181 if (quality == 0.0) | |
182 fscanf(sig, "%lf\n", &quality); | |
183 else | |
16
4987db85cfae
fix another scanf() warning
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
8
diff
changeset
|
184 fscanf(sig, "%*f\n"); |
0 | 185 if (threshold_noise == 0.0) |
186 fscanf(sig, "%lf\n", &threshold_noise); | |
187 else | |
16
4987db85cfae
fix another scanf() warning
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
8
diff
changeset
|
188 fscanf(sig, "%*f\n"); |
0 | 189 if (threshold_slope == 0.0) |
190 fscanf(sig, "%lf\n", &threshold_slope); | |
191 else | |
16
4987db85cfae
fix another scanf() warning
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
8
diff
changeset
|
192 fscanf(sig, "%*f\n"); |
0 | 193 if (blocksize == 0) |
194 fscanf(sig, "%d\n", &blocksize); | |
195 else | |
196 fscanf(sig, "%*d\n"); | |
197 fscanf(sig, "%d\n", &seed); | |
198 srandom(seed); | |
199 n_signature = NBITSTOBYTES(nbit_signature); | |
200 fread(signature, sizeof(char), n_signature, sig); | |
201 init_signature_bits(); | |
202 fscanf(sig, "\n"); | |
203 } | |
204 else { | |
205 fprintf(stderr, "%s: invalid signature file %s\n", progname, signature_name); | |
206 exit(1); | |
207 } | |
208 fclose(sig); | |
209 } | |
210 else { | |
211 fprintf(stderr, "%s: signature file not specified, use -s file option\n", progname); | |
212 exit(1); | |
213 } | |
214 | |
215 if (pattern1 <= 0 || pattern2 <= 0 || pattern1 > NPATTERN || pattern2 > NPATTERN) { | |
8 | 216 fprintf(stderr, "%s: invalid pattern type specified\n", progname); |
0 | 217 exit(1); |
218 } | |
219 | |
220 // read dimensions of input image file | |
221 pgm_readpgminit(in, &cols, &rows, &maxval, &format); | |
222 | |
223 // see if we can extract all signature bits | |
224 // we want at least half of the blocks untouched | |
225 if (((rows / blocksize) * (cols / blocksize)) < nbit_signature / 2) { | |
226 fprintf(stderr, "%s: image not large enough to contain %d bits of signature\n", progname, nbit_signature); | |
227 exit(1); | |
228 } | |
229 n_block = blocksize * blocksize; | |
230 | |
231 // allocate structure to remember which blocks we already touched, | |
232 // allow plenty of room to skip over blocks | |
8 | 233 if ((coords = alloc_coords(nbit_signature * 16)) == NULL) { |
0 | 234 fprintf(stderr, "%s: unable to allocate memory\n", progname); |
235 exit(1); | |
236 } | |
237 | |
238 // read in input image file | |
239 image = pgm_allocarray(cols, rows); | |
240 for (row = 0; row < rows; row++) | |
241 pgm_readpgmrow(in, image[row], cols, maxval, format); | |
242 | |
243 fclose(in); | |
244 | |
245 row = 0; | |
246 | |
247 // allocate memory for one block | |
248 block = alloc_grays(blocksize, blocksize); | |
249 | |
250 // allocate memory for zone classification | |
251 zone = alloc_grays(blocksize, blocksize); | |
252 | |
253 // allocate memory for category classification | |
254 category1 = alloc_grays(blocksize, blocksize); | |
255 category2 = alloc_grays(blocksize, blocksize); | |
256 | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
16
diff
changeset
|
257 // set up category classification array according to |
0 | 258 // pattern type parameter |
259 for (i = 0; i < blocksize; i++) | |
260 for (j = 0; j < blocksize; j++) { | |
261 category1[j][i] = lookup_pattern(pattern1, i, j); | |
262 category2[j][i] = lookup_pattern(pattern2, i, j); | |
263 } | |
264 | |
265 // allocate memory for slope calculation | |
266 slope = malloc(sizeof(double) * n_block); | |
267 | |
268 fprintf(out, "BRWM\n"); | |
269 fprintf(out, "%d\n", nbit_signature); | |
270 | |
271 // extract all the signature bits, one by one | |
272 n = 0; | |
273 while (n < nbit_signature) { | |
274 int xb; | |
275 int yb; | |
276 int blocktype; | |
277 double smax; | |
278 int alpha, beta_minus, beta_plus; | |
279 double mean_1A, mean_1B, mean_2A, mean_2B, mean_1, mean_2; | |
280 int n_1A, n_1B, n_2A, n_2B, n_1, n_2; | |
281 double sigma, sigma_1, sigma_2; | |
282 int zone1_ok, zone2_ok; | |
283 | |
284 // find an unused block randomly, depending on seed | |
285 do { | |
286 xb = random() % (cols / blocksize); | |
287 yb = random() % (rows / blocksize); | |
288 } while (add_coord(coords, xb, yb) < 0); | |
289 | |
290 // copy image block | |
8 | 291 fprintf(stderr, "XXX1 %d %d %d\n", xb*blocksize, yb*blocksize, blocksize); |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
16
diff
changeset
|
292 copy_grays_to_block(block, image, xb * blocksize, yb * blocksize, blocksize, blocksize); |
8 | 293 fprintf(stderr, "XXX2\n"); |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
16
diff
changeset
|
294 |
0 | 295 if (verbose > 0) |
296 fprintf(stderr, "detecting bit #%d in block at (%d/%d)\n", n, xb * blocksize, yb * blocksize); | |
297 | |
298 // sort luminance values in block to represent increasing function F | |
299 sort_grays(block[0], n_block); | |
300 | |
301 // calculate slopes of F and determine smax, the max. slope of F | |
302 // the index where smax occures is called alpha | |
303 alpha = 0; | |
304 smax = 0.0; | |
305 for (i = 0; i < n_block - 1; i++) { | |
306 slope[i] = block[0][i + 1] - block[0][i]; | |
307 if (slope[i] > smax) { | |
308 smax = slope[i]; | |
309 alpha = i; | |
310 } | |
311 } | |
312 slope[n_block - 1] = 0; | |
313 | |
314 // block type classification | |
315 blocktype = BLOCKTYPE_UNKNOWN; | |
316 | |
317 if (smax < threshold_noise) { | |
318 // block has noise contrast | |
319 beta_minus = beta_plus = alpha; | |
320 blocktype = BLOCKTYPE_NOISE; | |
321 } | |
322 else { | |
323 // block has progressive or hard contrast, let's find out... | |
324 | |
325 beta_minus = alpha - 1; | |
326 while (beta_minus >= 0 && smax - slope[beta_minus] <= threshold_slope) | |
327 beta_minus--; | |
328 | |
329 beta_plus = alpha + 1; | |
330 while (beta_plus < n_block && smax - slope[beta_plus] <= threshold_slope) | |
331 beta_plus++; | |
332 | |
333 if (beta_minus + 1 == alpha && beta_plus - 1 == alpha) | |
334 blocktype = BLOCKTYPE_HARD; | |
335 else | |
336 blocktype = BLOCKTYPE_PROGRESSIVE; | |
337 } | |
338 | |
339 if (verbose > 1) { | |
340 fprintf(stderr, "blocktype: %d\n", blocktype); | |
341 fprintf(stderr, "Smax = %lf, alpha = %d, beta- = %d, beta+ = %d\n", smax, alpha, beta_minus, beta_plus); | |
342 } | |
343 | |
344 // block pixel classification | |
345 for (i = 0; i < blocksize; i++) | |
346 for (j = 0; j < blocksize; j++) { | |
347 gray pixel = image[yb * blocksize + j][xb * blocksize + i]; | |
348 zone[j][i] = ZONE_VOID; | |
349 switch (blocktype) { | |
350 case BLOCKTYPE_PROGRESSIVE: | |
351 case BLOCKTYPE_HARD: | |
352 if (pixel < block[0][beta_minus]) | |
353 zone[j][i] = ZONE_1; | |
354 else if (pixel > block[0][beta_plus]) | |
355 zone[j][i] = ZONE_2; | |
356 break; | |
357 case BLOCKTYPE_NOISE: | |
358 if (pixel < block[0][n_block / 2]) | |
359 zone[j][i] = ZONE_1; | |
360 else if (pixel > block[0][n_block / 2]) | |
361 zone[j][i] = ZONE_2; | |
362 break; | |
363 default: | |
364 fprintf(stderr, "%s: invalid block type\n", progname); | |
365 break; | |
366 } | |
367 } | |
368 | |
369 // calculate mean values for zone/categories | |
370 mean_1A = mean_1B = mean_2A = mean_2B = mean_1 = mean_2 = 0.0; | |
371 n_1A = n_1B = n_2A = n_2B = n_1 = n_2 = 0; | |
372 for (i = 0; i < blocksize; i++) | |
373 for (j = 0; j < blocksize; j++) { | |
374 gray pixel = image[yb * blocksize + j][xb * blocksize + i]; | |
375 int pixel_zone = zone[j][i]; | |
376 int pixel_category = CATEGORY_VOID; | |
377 if (pixel_zone == ZONE_1) | |
378 pixel_category = category1[j][i]; | |
379 else if (pixel_zone == ZONE_2) | |
380 pixel_category = category2[j][i]; | |
381 | |
382 switch (pixel_zone | pixel_category) { | |
383 case CLASSIFICATION_1A: | |
384 n_1++; | |
385 n_1A++; | |
386 mean_1A += pixel; | |
387 mean_1 += pixel; | |
388 break; | |
389 case CLASSIFICATION_1B: | |
390 n_1++; | |
391 n_1B++; | |
392 mean_1B += pixel; | |
393 mean_1 += pixel; | |
394 break; | |
395 case CLASSIFICATION_2A: | |
396 n_2++; | |
397 n_2A++; | |
398 mean_2A += pixel; | |
399 mean_2 += pixel; | |
400 break; | |
401 case CLASSIFICATION_2B: | |
402 n_2++; | |
403 n_2B++; | |
404 mean_2B += pixel; | |
405 mean_2 += pixel; | |
406 break; | |
407 } | |
408 } | |
409 | |
410 if (n_1 && n_1A && n_1B) { | |
411 mean_1 /= (double) n_1; | |
412 mean_1A /= (double) n_1A; | |
413 mean_1B /= (double) n_1B; | |
414 zone1_ok = 1; | |
415 } | |
416 else { | |
417 mean_1 = mean_1A = mean_1B = 0.0; | |
418 zone1_ok = 0; | |
419 if (verbose > 0) | |
420 fprintf(stderr, "zone 1 unusable\n"); | |
421 } | |
422 | |
423 if (n_2 && n_2A && n_2B) { | |
424 mean_2 /= (double) n_2; | |
425 mean_2A /= (double) n_2A; | |
426 mean_2B /= (double) n_2B; | |
427 zone2_ok = 1; | |
428 } | |
429 else { | |
430 mean_2 = mean_2A = mean_2B = 0.0; | |
431 zone2_ok = 0; | |
432 if (verbose > 0) | |
433 fprintf(stderr, "zone 2 unusable\n"); | |
434 } | |
435 | |
436 // bit extraction | |
437 if (zone1_ok && zone2_ok) { | |
438 sigma_1 = mean_1A - mean_1B; | |
439 sigma_2 = mean_2A - mean_2B; | |
440 | |
441 if (verbose > 2) { | |
442 fprintf(stderr, "m_1A = %lf, m_1B = %lf\n", mean_1A, mean_1B); | |
443 fprintf(stderr, "m_2A = %lf, m_2B = %lf\n", mean_2A, mean_2B); | |
444 fprintf(stderr, "sigma1 = %lf, sigma2 = %lf\n", sigma_1, sigma_2); | |
445 } | |
446 | |
447 #define EPSILON 0.001 | |
448 if (fabs(sigma_1 * sigma_2) < EPSILON) { | |
449 // case 3 | |
450 sigma = MAX(fabs(sigma_1), fabs(sigma_2)); | |
451 set_signature_bit(n, sigma > 0.0); | |
452 if (verbose > 0) | |
453 fprintf(stderr, "case 3, bit #%d = %d\n", n, sigma > 0.0); | |
454 } | |
455 else if (sigma_1 * sigma_2 > 0.0) { | |
456 // case 1 | |
457 set_signature_bit(n, sigma_1 > 0.0); | |
458 if (verbose > 0) | |
459 fprintf(stderr, "case 1, bit #%d = %d\n", n, sigma_1 > 0.0); | |
460 } | |
461 else if (sigma_1 * sigma_2 < 0.0) { | |
462 // case 2 | |
463 sigma = (double) (n_1A + n_1B) * sigma_1 + (double) (n_2A + n_2B) * sigma_2; | |
464 set_signature_bit(n, sigma > 0.0); | |
465 if (verbose > 0) | |
466 fprintf(stderr, "case 2, bit #%d = %d\n", n, sigma > 0.0); | |
467 } | |
468 } | |
469 else if (zone1_ok) { | |
470 set_signature_bit(n, mean_1A > mean_1B); | |
471 if (verbose > 0) | |
472 fprintf(stderr, "case 4, bit #%d = %d\n", n, mean_1A > mean_1B); | |
473 } | |
474 else if (zone2_ok) { | |
475 set_signature_bit(n, mean_2A > mean_2B); | |
476 if (verbose > 0) | |
477 fprintf(stderr, "case 5, bit #%d = %d\n", n, mean_2A > mean_2B); | |
478 } | |
479 else { | |
480 // pathological case - can it ever happen? | |
481 if (verbose > 0) | |
482 fprintf(stderr, "block skipped\n"); | |
483 if (!skipping) continue; | |
484 } | |
485 | |
486 n++; | |
487 } | |
488 | |
489 free_grays(category2); | |
490 free_grays(category1); | |
491 free_grays(zone); | |
492 free_grays(block); | |
493 | |
494 // write extracted signature | |
495 | |
496 fwrite(signature, sizeof(char), n_signature, out); | |
497 fclose(out); | |
498 | |
499 pgm_freearray(image, rows); | |
500 | |
501 exit(0); | |
502 } |