Mercurial > hg > wm
comparison Meerwald/wm_koch_e.c @ 18:3bdb67e76858
mse opt.
author | Peter Meerwald <pmeerw@cosy.sbg.ac.at> |
---|---|
date | Fri, 30 Jan 2009 12:46:49 +0100 |
parents | 824d192e5614 |
children |
comparison
equal
deleted
inserted
replaced
17:824d192e5614 | 18:3bdb67e76858 |
---|---|
1 #include "wm.h" | 1 #include "wm.h" |
2 #include "dct.h" | 2 #include "dct.h" |
3 #include "signature.h" | 3 #include "signature.h" |
4 #include "coord.h" | 4 #include "coord.h" |
5 #include "gray.h" | |
5 #include "pgm.h" | 6 #include "pgm.h" |
6 | 7 |
7 char *progname; | 8 char *progname; |
9 | |
10 double sign(double x) { | |
11 if (x >= 0.0) return 1.0; | |
12 else return -1.0; | |
13 } | |
14 | |
15 double try_modif(gray **image_block, double **dcts, int c1, int c2, double w1, double w2) { | |
16 int i, j; | |
17 gray **altered_block; | |
18 double **altered_dcts; | |
19 double sum; | |
20 | |
21 altered_block = alloc_grays_8x8(); | |
22 altered_dcts = alloc_coeffs_8x8(); | |
23 | |
24 for (i = 0; i < 8; i++) { | |
25 memcpy(altered_dcts[i], dcts[i], sizeof(double) * 8); | |
26 } | |
27 | |
28 // put the changed coefficients back to black | |
29 altered_dcts[c1 / NJPEG][c1 % NJPEG] = w1; | |
30 altered_dcts[c2 / NJPEG][c2 % NJPEG] = w2; | |
31 | |
32 dequantize_8x8(altered_dcts); | |
33 | |
34 idct_block_8x8(altered_dcts, altered_block, 0, 0); | |
35 | |
36 // compute MSE | |
37 sum = 0.0; | |
38 for (i = 0; i < 8; i++) { | |
39 for (j = 0; j < 8; j++) { | |
40 double ib = image_block[i][j]; | |
41 double ab = altered_block[i][j]; | |
42 sum += (ib - ab) * (ib - ab); | |
43 } | |
44 } | |
45 sum /= 64.0; | |
46 | |
47 free(altered_block); | |
48 free(altered_dcts); | |
49 | |
50 return sum; | |
51 } | |
8 | 52 |
9 void usage(void) { | 53 void usage(void) { |
10 fprintf(stderr, "usage: %s [-h] [-l n] [-o file] [-q n] [-v n] -s file file\n", progname); | 54 fprintf(stderr, "usage: %s [-h] [-l n] [-o file] [-q n] [-v n] -s file file\n", progname); |
11 fprintf(stderr, "\t-h\t\tprint usage\n"); | 55 fprintf(stderr, "\t-h\t\tprint usage\n"); |
12 fprintf(stderr, "\t-l n\t\tsignature robustness factor\n"); | 56 fprintf(stderr, "\t-l n\t\tsignature robustness factor\n"); |
42 | 86 |
43 struct coords *coords; | 87 struct coords *coords; |
44 | 88 |
45 gray **image; | 89 gray **image; |
46 double **dcts; | 90 double **dcts; |
91 gray **image_block; | |
47 | 92 |
48 progname = argv[0]; | 93 progname = argv[0]; |
49 | 94 |
50 pgm_init(&argc, argv); wm_init(); | 95 pgm_init(&argc, argv); wm_init(); |
51 | 96 |
159 | 204 |
160 init_dct_8x8(); | 205 init_dct_8x8(); |
161 init_quantum_JPEG_lumin(quantization); | 206 init_quantum_JPEG_lumin(quantization); |
162 | 207 |
163 dcts = alloc_coeffs_8x8(); | 208 dcts = alloc_coeffs_8x8(); |
209 image_block = alloc_grays_8x8(); | |
164 | 210 |
165 if ((coords = alloc_coords(nbit_signature)) == NULL) { | 211 if ((coords = alloc_coords(nbit_signature)) == NULL) { |
166 fprintf(stderr, "%s: unable to allocate memory\n", progname); | 212 fprintf(stderr, "%s: unable to allocate memory\n", progname); |
167 exit(1); | 213 exit(1); |
168 } | 214 } |
181 int xb; | 227 int xb; |
182 int yb; | 228 int yb; |
183 int c1, c2; | 229 int c1, c2; |
184 double v1, v2; | 230 double v1, v2; |
185 double w1, w2; | 231 double w1, w2; |
186 double diff, abs_diff; | 232 double best_w1, best_w2; |
233 double diff; | |
234 double mod; | |
235 double try; | |
236 double best_mse; | |
237 int no_mse_opt = 0; | |
187 | 238 |
188 // randomly select a block, check to get distinct blocks | 239 // randomly select a block, check to get distinct blocks |
189 // (don't watermark a block twice) | 240 // (don't watermark a block twice) |
190 do { | 241 do { |
191 xb = random() % (cols / NJPEG); | 242 xb = random() % (cols / NJPEG); |
193 } while (add_coord(coords, xb, yb) < 0); | 244 } while (add_coord(coords, xb, yb) < 0); |
194 | 245 |
195 // do the forward 8x8 DCT of that block | 246 // do the forward 8x8 DCT of that block |
196 fdct_block_8x8(image, xb * NJPEG, yb * NJPEG, dcts); | 247 fdct_block_8x8(image, xb * NJPEG, yb * NJPEG, dcts); |
197 | 248 |
249 copy_grays_to_block(image_block, image, xb*NJPEG, yb*NJPEG, NJPEG, NJPEG); | |
250 | |
198 // randomly select two distinct coefficients from block | 251 // randomly select two distinct coefficients from block |
199 // only accept coefficients in the middle frequency range | 252 // only accept coefficients in the middle frequency range |
200 do { | 253 do { |
201 c1 = (random() % (NJPEG * NJPEG - 2)) + 1; | 254 c1 = (random() % (NJPEG * NJPEG - 2)) + 1; |
202 c2 = (random() % (NJPEG * NJPEG - 2)) + 1; | 255 c2 = (random() % (NJPEG * NJPEG - 2)) + 1; |
211 print_coeffs_8x8(dcts); | 264 print_coeffs_8x8(dcts); |
212 | 265 |
213 v1 = dcts[c1 / NJPEG][c1 % NJPEG]; | 266 v1 = dcts[c1 / NJPEG][c1 % NJPEG]; |
214 v2 = dcts[c2 / NJPEG][c2 % NJPEG]; | 267 v2 = dcts[c2 / NJPEG][c2 % NJPEG]; |
215 | 268 |
269 best_w1 = DBL_MAX, best_w2 = DBL_MAX; | |
270 try = 0.0; | |
271 best_mse = DBL_MAX; | |
272 | |
216 diff = fabs(v1) - fabs(v2); | 273 diff = fabs(v1) - fabs(v2); |
217 abs_diff = (fabs(diff) + quality) / 2.0; | 274 |
218 | 275 if (get_signature_bit(n)) |
219 // modify coefficient's relationship to embed signature bit | 276 mod = fabs(quality - ( fabs(v1) - fabs(v2) )); |
220 // using mean square error to minimize error | 277 else |
221 if (get_signature_bit(n)) { | 278 mod = fabs(quality - (fabs(v2) - fabs(v1))); |
222 if (diff < quality) { | 279 |
223 // we have to impose the relationship, does not occur naturally | 280 if (verbose > 2) |
224 w1 = (v1 > 0.0) ? (v1 + abs_diff) : (v1 - abs_diff); | 281 fprintf(stderr, "%d / %d: %.2f %.2f %.2f | %d\n", xb, yb, diff, v1, v2, get_signature_bit(n)); |
225 w2 = v2; | 282 |
226 } | 283 while (try <= mod) { |
227 else { | |
228 w1 = v1; | 284 w1 = v1; |
229 w2 = v2; | 285 w2 = v2; |
230 } | 286 |
231 } | 287 // modify coefficient's relationship to embed signature bit |
232 else { | 288 // using mean square error to minimize error |
233 if (diff > -quality) { | 289 if (get_signature_bit(n)) { |
234 // force the relationship | 290 if (diff < quality) { |
235 w1 = v1; | 291 // we have to impose the relationship, does not occur naturally |
236 w2 = (v2 > 0.0) ? (v2 + abs_diff) : (v2 - abs_diff); | 292 w1 = sign(v1)*(fabs(v1) + mod - try); |
237 } | 293 w2 = sign(v2)*(fabs(v2) - try); |
238 else { | 294 } |
239 w1 = v1; | 295 } |
240 w2 = v2; | 296 else { |
241 } | 297 if (diff > -quality) { |
298 // force the relationship | |
299 w2 = sign(v2)*(fabs(v2) + mod - try); | |
300 w1 = sign(v1)*(fabs(v1) - try); | |
301 } | |
302 } | |
303 | |
304 double mse = try_modif(image_block, dcts, c1, c2, w1, w2); | |
305 if (mse < best_mse) { | |
306 best_w1 = w1; | |
307 best_w2 = w2; | |
308 best_mse = mse; | |
309 } | |
310 | |
311 if (verbose > 2) | |
312 fprintf(stderr, "%d / %d: MSE %.2f %.2f; %.2f: %.2f %.2f\n", xb, yb, mse, best_mse, try, w1, w2); | |
313 | |
314 if (fabs(mse) == 1e-3) | |
315 break; | |
316 | |
317 if (fabs(fabs(w1) - fabs(w2) + quality) > 1e-3) | |
318 break; | |
319 | |
320 if (no_mse_opt) | |
321 break; | |
322 | |
323 try += 0.05; | |
242 } | 324 } |
243 | 325 |
244 if (verbose > 1) | 326 if (verbose > 1) |
245 fprintf(stderr, " %f -> %f, %f -> %f\n", v1, w1, v2, w2); | 327 fprintf(stderr, " %f -> %f, %f -> %f\n", v1, best_w1, v2, best_w2); |
246 | 328 |
247 // put the changed coefficients back to black | 329 // put the changed coefficients back to black |
248 dcts[c1 / NJPEG][c1 % NJPEG] = w1; | 330 dcts[c1 / NJPEG][c1 % NJPEG] = best_w1; |
249 dcts[c2 / NJPEG][c2 % NJPEG] = w2; | 331 dcts[c2 / NJPEG][c2 % NJPEG] = best_w2; |
250 | 332 |
251 // the obvious :-) | 333 // the obvious :-) |
252 dequantize_8x8(dcts); | 334 dequantize_8x8(dcts); |
253 | 335 |
254 // do the inverse DCT on the modified 8x8 block | 336 // do the inverse DCT on the modified 8x8 block |
255 idct_block_8x8(dcts, image, xb * NJPEG, yb * NJPEG); | 337 idct_block_8x8(dcts, image, xb * NJPEG, yb * NJPEG); |
256 | 338 |
257 n++; | 339 n++; |
258 } | 340 } |
259 | 341 |
342 free_grays(image_block); | |
260 free_coeffs(dcts); | 343 free_coeffs(dcts); |
261 | 344 |
262 pgm_writepgminit(out, cols, rows, maxval, 0); | 345 pgm_writepgminit(out, cols, rows, maxval, 0); |
263 for (row = 0; row < rows; row++) | 346 for (row = 0; row < rows; row++) |
264 pgm_writepgmrow(out, image[row], cols, maxval, 0); | 347 pgm_writepgmrow(out, image[row], cols, maxval, 0); |