Mercurial > hg > wm
annotate Meerwald/wm_bruyn_d.c @ 22:d8551fb39a5e default tip
Added tag v0.6 for changeset 1c4ccd635a68
| author | Peter Meerwald-Stadler <pmeerw@pmeerw.net> |
|---|---|
| date | Sat, 28 Jan 2023 23:57:51 +0100 |
| parents | bd669312f068 |
| 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" | |
| 7 #include "pgm.h" | |
| 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 } |