--- a/Meerwald/cmp_koch_sig.c	Thu Jan 29 18:28:05 2009 +0100
+++ b/Meerwald/cmp_koch_sig.c	Fri Jan 30 12:46:49 2009 +0100
@@ -132,6 +132,15 @@
     else
       corr--;
 
+  if (verbose > 1) {
+    for (i = 0; i < nbit_signature2; i++)
+      fprintf(stderr, "%d", get_signature1_bit(i % nbit_signature1));
+    fprintf(stderr, "\n");
+    for (i = 0; i < nbit_signature2; i++)
+      fprintf(stderr, "%d", get_signature2_bit(i));
+    fprintf(stderr, "\n");
+  }
+
   if (correlation_only)
     fprintf(out, "%lf\n", (double) corr / nbit_signature2);
   else {

--- a/Meerwald/wm_koch_e.c	Thu Jan 29 18:28:05 2009 +0100
+++ b/Meerwald/wm_koch_e.c	Fri Jan 30 12:46:49 2009 +0100
@@ -2,10 +2,54 @@
 #include "dct.h"
 #include "signature.h"
 #include "coord.h"
+#include "gray.h"
 #include "pgm.h"
 
 char *progname;
 
+double sign(double x) {
+    if (x >= 0.0) return 1.0;
+    else return -1.0;
+}
+
+double try_modif(gray **image_block, double **dcts, int c1, int c2, double w1, double w2) {
+    int i, j;
+    gray **altered_block;
+    double **altered_dcts;
+    double sum;
+
+    altered_block = alloc_grays_8x8();
+    altered_dcts = alloc_coeffs_8x8();
+
+    for (i = 0; i < 8; i++) {
+        memcpy(altered_dcts[i], dcts[i], sizeof(double) * 8);
+    }
+        
+    // put the changed coefficients back to black
+    altered_dcts[c1 / NJPEG][c1 % NJPEG] = w1;
+    altered_dcts[c2 / NJPEG][c2 % NJPEG] = w2;
+
+    dequantize_8x8(altered_dcts);
+
+    idct_block_8x8(altered_dcts, altered_block, 0, 0);
+    
+    // compute MSE
+    sum = 0.0;
+    for (i = 0; i < 8; i++) {
+        for (j = 0; j < 8; j++) {  
+            double ib = image_block[i][j];
+            double ab = altered_block[i][j];
+            sum += (ib - ab) * (ib - ab);
+        }
+    }
+    sum /= 64.0;
+    
+    free(altered_block);
+    free(altered_dcts);
+    
+    return sum;
+}
+
 void usage(void) {
   fprintf(stderr, "usage: %s [-h] [-l n] [-o file] [-q n] [-v n] -s file file\n", progname);
   fprintf(stderr, "\t-h\t\tprint usage\n");
@@ -44,6 +88,7 @@
 
   gray **image;
   double **dcts;
+  gray **image_block;
 
   progname = argv[0]; 
 
@@ -161,6 +206,7 @@
   init_quantum_JPEG_lumin(quantization);
 
   dcts = alloc_coeffs_8x8();
+  image_block = alloc_grays_8x8();
 
   if ((coords = alloc_coords(nbit_signature)) == NULL) {
     fprintf(stderr, "%s: unable to allocate memory\n", progname);
@@ -183,7 +229,12 @@
     int c1, c2;
     double v1, v2;
     double w1, w2;
-    double diff, abs_diff;
+    double best_w1, best_w2;
+    double diff;
+    double mod;
+    double try;
+    double best_mse;
+    int no_mse_opt = 0;
 
     // randomly select a block, check to get distinct blocks
     // (don't watermark a block twice)
@@ -195,6 +246,8 @@
     // do the forward 8x8 DCT of that block
     fdct_block_8x8(image, xb * NJPEG, yb * NJPEG, dcts);
 
+    copy_grays_to_block(image_block, image, xb*NJPEG, yb*NJPEG, NJPEG, NJPEG);
+
     // randomly select two distinct coefficients from block
     // only accept coefficients in the middle frequency range
     do {
@@ -213,40 +266,69 @@
     v1 = dcts[c1 / NJPEG][c1 % NJPEG];
     v2 = dcts[c2 / NJPEG][c2 % NJPEG];
 
+    best_w1 = DBL_MAX, best_w2 = DBL_MAX;
+    try = 0.0;
+    best_mse = DBL_MAX;
+
     diff = fabs(v1) - fabs(v2);
-    abs_diff = (fabs(diff) + quality) / 2.0;
 
-    // modify coefficient's relationship to embed signature bit
-    // using mean square error to minimize error
-    if (get_signature_bit(n)) {
-      if (diff < quality) {
-        // we have to impose the relationship, does not occur naturally
-        w1 = (v1 > 0.0) ? (v1 + abs_diff) : (v1 - abs_diff);
-        w2 = v2;
-      }
-      else {
+    if (get_signature_bit(n)) 
+      mod = fabs(quality - ( fabs(v1) - fabs(v2) ));
+    else
+      mod = fabs(quality -  (fabs(v2) - fabs(v1)));
+    
+    if (verbose > 2)
+        fprintf(stderr, "%d / %d: %.2f %.2f %.2f | %d\n", xb, yb, diff, v1, v2, get_signature_bit(n));    
+    
+    while (try <= mod) {
         w1 = v1;
         w2 = v2;
-      }
-    }
-    else {
-      if (diff > -quality) {
-        // force the relationship
-        w1 = v1;
-        w2 = (v2 > 0.0) ? (v2 + abs_diff) : (v2 - abs_diff);
-      }
-      else {
-        w1 = v1;
-        w2 = v2;
-      }
+
+        // modify coefficient's relationship to embed signature bit
+        // using mean square error to minimize error
+        if (get_signature_bit(n)) {
+          if (diff < quality) {
+            // we have to impose the relationship, does not occur naturally
+            w1 = sign(v1)*(fabs(v1) + mod - try);
+            w2 = sign(v2)*(fabs(v2) - try);
+          }
+        }
+        else {
+          if (diff > -quality) {
+            // force the relationship
+            w2 = sign(v2)*(fabs(v2) + mod - try);
+            w1 = sign(v1)*(fabs(v1) - try);
+          }
+        }
+
+        double mse = try_modif(image_block, dcts, c1, c2, w1, w2);
+        if (mse < best_mse) {
+            best_w1 = w1;
+            best_w2 = w2;
+            best_mse = mse;
+        }
+        
+        if (verbose > 2)
+            fprintf(stderr, "%d / %d: MSE %.2f %.2f; %.2f: %.2f %.2f\n", xb, yb, mse, best_mse, try, w1, w2);
+        
+        if (fabs(mse) == 1e-3)
+            break;
+        
+        if (fabs(fabs(w1) - fabs(w2) + quality) > 1e-3)
+            break;
+
+        if (no_mse_opt)
+            break;
+
+        try += 0.05;
     }
 
     if (verbose > 1)
-      fprintf(stderr, "  %f -> %f, %f -> %f\n", v1, w1, v2, w2);
+      fprintf(stderr, "  %f -> %f, %f -> %f\n", v1, best_w1, v2, best_w2);
 
     // put the changed coefficients back to black
-    dcts[c1 / NJPEG][c1 % NJPEG] = w1;
-    dcts[c2 / NJPEG][c2 % NJPEG] = w2;
+    dcts[c1 / NJPEG][c1 % NJPEG] = best_w1;
+    dcts[c2 / NJPEG][c2 % NJPEG] = best_w2;
 
     // the obvious :-)
     dequantize_8x8(dcts);
@@ -257,6 +339,7 @@
     n++;
   }
 
+  free_grays(image_block);
   free_coeffs(dcts);
 
   pgm_writepgminit(out, cols, rows, maxval, 0);