Mercurial > hg > wm

comparison Meerwald/wm_koch_e.c @ 18:3bdb67e76858

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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);

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.