stc/exp/ml_stc_linux_make_v1.0/ml_stc_src/coding_loss.cpp

   1 #include <iostream>\r
   2 #include <iomanip>\r
   3 #include <cmath>\r
   4 #include "coding_loss.h"\r
   5 \r
   6 // binary entropy function in double precission\r
   7 inline double binary_entropy(double x) {\r
   8 \r
   9     double const LOG2 = log(2.0);\r
  10     double const EPS = std::numeric_limits<double>::epsilon();\r
  11     double z;\r
  12 \r
  13     if ((x<EPS) || ((1-x)<EPS)) {\r
  14         return 0;\r
  15     } else {\r
  16         z = (-x*log(x)-(1-x)*log(1-x))/LOG2;\r
  17         return z;\r
  18     }\r
  19 }\r
  20 \r
  21 double calc_coding_loss(const u8 *cover, unsigned int cover_length, unsigned int message_length, const double *costs, const u8 *stego) {\r
  22     \r
  23     unsigned int i;\r
  24     double avg_distortion = 0; // average distortion between cover and stego\r
  25     double opt_rel_payload = 0;\r
  26     double beta, actual_rel_payload, coding_loss;\r
  27 \r
  28     // calculate average distortion per cover element\r
  29     for(i=0;i<cover_length;i++) { avg_distortion += (cover[i]==stego[i]) ? 0 : costs[i]; }\r
  30     avg_distortion /= cover_length;\r
  31 \r
  32     // find optimal value of parameter beta such that optimal coding achieves avg_distortion\r
  33     beta = calculate_beta_from_distortion(costs, cover_length, avg_distortion);\r
  34     // calculate optimal relative payload achieved by optimal coding method\r
  35     for(i=0;i<cover_length;i++) {\r
  36         double x = exp(-beta*costs[i]);\r
  37         opt_rel_payload += binary_entropy(x/(1.0+x)); // measure entropy of the optimal flipping noise\r
  38     }\r
  39     opt_rel_payload /= cover_length;\r
  40     actual_rel_payload = ((double)message_length)/((double)cover_length);\r
  41 \r
  42     coding_loss = 1.0 - actual_rel_payload/opt_rel_payload;\r
  43     return coding_loss;\r
  44 }\r
  45 \r
  46 /*\r
  47   find beta such that average distortion caused by optimal flipping noise is avg_distortion\r
  48   using binary search\r
  49 */\r
  50 double calculate_beta_from_distortion(const double *costs, unsigned int n, double avg_distortion) {\r
  51 \r
  52     double dist1, dist2, dist3, beta1, beta2, beta3;\r
  53     int j = 0;\r
  54     unsigned int i;\r
  55     const double INF = std::numeric_limits<double>::infinity();\r
  56 \r
  57     beta2 = -1; // initial value - error\r
  58     beta1 = 0; dist1 = 0;\r
  59     for(i=0;i<n;i++) {\r
  60         dist1 += (costs[i]==INF) ? 0 : costs[i];\r
  61     }\r
  62     dist1 /= n*2;\r
  63     beta3 = 1e+5; dist3 = avg_distortion+1; // this is just an initial value\r
  64     while(dist3>avg_distortion) {\r
  65         beta3 *= 2;\r
  66         for(i=0;i<n;i++) {\r
  67             double p_flip = 1.0-1.0/(1.0+exp(-beta3*costs[i]));\r
  68             if (costs[i]!=INF) { dist3 += costs[i]*p_flip; }\r
  69         }\r
  70         dist3 /= n;\r
  71         j++;\r
  72         if (j>100) {\r
  73             // beta is probably unbounded => it seems that we cannot find beta such that\r
  74             // relative distortion will be smaller than requested. Binary search does not make sense here.\r
  75             return -1;\r
  76         }\r
  77     }\r
  78     while ( (dist1-dist3>(avg_distortion/1000.0)) ) { // iterative search for parameter beta\r
  79         beta2 = beta1+(beta3-beta1)/2; dist2 = 0;\r
  80         for(i=0;i<n;i++) {\r
  81             if (costs[i]!=INF) {\r
  82                 double p_flip = 1.0-1.0/(1.0+exp(-beta2*costs[i]));\r
  83                 if (costs[i]!=INF) { dist2 += costs[i]*p_flip; }\r
  84             }\r
  85         }\r
  86         dist2 /= n;\r
  87         if (dist2<avg_distortion) {\r
  88             beta3 = beta2;\r
  89             dist3 = dist2;\r
  90         } else {\r
  91             beta1 = beta2;\r
  92             dist1 = dist2;\r
  93         }\r
  94     }\r
  95 \r
  96     return beta2;\r
  97 }\r