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

   1 //\r
   2 // MEX-File: dv_embed\r
   3 // Usage:\r
   4 // [double scalar price, uint8 vector stego] = dv_embed(uint8 vector cover, uint8 vector message, uint8/double vector pricevec, int scalar constrLength = 10)\r
   5 \r
   6 typedef unsigned int u32;\r
   7 typedef unsigned short u16;\r
   8 typedef unsigned char u8;\r
   9 \r
  10 #include <cstdlib>\r
  11 #include <cstring>\r
  12 #include <cmath>\r
  13 #include <cfloat>\r
  14 #include <limits>\r
  15 #include <emmintrin.h>\r
  16 #include "mex.h"\r
  17 #include "common.h"\r
  18 \r
  19 #define MATLAB_in_cover 0\r
  20 #define MATLAB_in_message 1\r
  21 #define MATLAB_in_pricevec 2\r
  22 #define MATLAB_in_constrLength 3\r
  23 \r
  24 #define MATLAB_out_price 0\r
  25 #define MATLAB_out_stego 1\r
  26 \r
  27 void *aligned_malloc(unsigned int bytes, int align)\r
  28 {\r
  29         int shift;\r
  30         char *temp = (char *)malloc(bytes + align);\r
  31 \r
  32         if(temp == NULL)\r
  33                 return temp;\r
  34         shift = align - (int)(((unsigned long long)temp) & (align - 1));\r
  35         temp = temp + shift;\r
  36         temp[-1] = shift;\r
  37         return (void *)temp;\r
  38 }\r
  39 \r
  40 void aligned_free(void *vptr)\r
  41 {\r
  42         char *ptr = (char *)vptr;\r
  43         free(ptr - ptr[-1]);\r
  44         return;\r
  45 }\r
  46 \r
  47 inline __m128i maxLessThan255(const __m128i v1, const __m128i v2)\r
  48 {\r
  49         register __m128i mask = _mm_set1_epi32(0xffffffff);\r
  50         return _mm_max_epu8(_mm_andnot_si128(_mm_cmpeq_epi8(v1, mask), v1), _mm_andnot_si128(_mm_cmpeq_epi8(v2, mask), v2));\r
  51 }\r
  52 \r
  53 inline u8 max16B(__m128i maxp)\r
  54 {\r
  55         u8 mtemp[4];\r
  56         maxp = _mm_max_epu8(maxp, _mm_srli_si128(maxp, 8));\r
  57         maxp = _mm_max_epu8(maxp, _mm_srli_si128(maxp, 4));\r
  58         *((int*)mtemp) = _mm_cvtsi128_si32(maxp);\r
  59         if(mtemp[2] > mtemp[0])\r
  60                 mtemp[0] = mtemp[2];\r
  61         if(mtemp[3] > mtemp[1])\r
  62                 mtemp[1] = mtemp[3];\r
  63         if(mtemp[1] > mtemp[0])\r
  64                 return mtemp[1];\r
  65         else \r
  66                 return mtemp[0];\r
  67 }\r
  68 \r
  69 inline u8 min16B(__m128i minp)\r
  70 {\r
  71         u8 mtemp[4];\r
  72         minp = _mm_min_epu8(minp, _mm_srli_si128(minp, 8));\r
  73         minp = _mm_min_epu8(minp, _mm_srli_si128(minp, 4));\r
  74         *((int*)mtemp) = _mm_cvtsi128_si32(minp);\r
  75         if(mtemp[2] < mtemp[0])\r
  76                 mtemp[0] = mtemp[2];\r
  77         if(mtemp[3] < mtemp[1])\r
  78                 mtemp[1] = mtemp[3];\r
  79         if(mtemp[1] < mtemp[0])\r
  80                 return mtemp[1];\r
  81         else \r
  82                 return mtemp[0];\r
  83 }\r
  84 \r
  85 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
  86 {\r
  87         int matrixheight, vectorlength, height, i, j, k, l, index, index2, parts, m, sseheight, altm, pathindex, syndromelength, matrixwidth;\r
  88         u32 column, colmask, state;\r
  89         double totalprice;\r
  90         bool usefloat;\r
  91 \r
  92         void *temp, *pricevectorv;\r
  93         u8 *vector, *syndrome, *closestData, *ssedone, *stepdown;\r
  94         u32 *path, *columns[2];\r
  95         int *matrices, *widths;\r
  96 \r
  97         if(nrhs < 3) {\r
  98                 mexErrMsgTxt("Too few input parameters (3 or 4 expected)");\r
  99                 return;\r
 100         }\r
 101         \r
 102         if(!mxIsUint8(prhs[MATLAB_in_cover])) {\r
 103                 mexErrMsgTxt("The input vector must be of type uint8.");\r
 104                 return;\r
 105         }\r
 106         vector = (u8*)mxGetPr(prhs[MATLAB_in_cover]);\r
 107         vectorlength = (int)mxGetM(prhs[MATLAB_in_cover]);\r
 108         \r
 109         if(!mxIsUint8(prhs[MATLAB_in_message])) {\r
 110                 mexErrMsgTxt("The syndrome must be of type uint8.");\r
 111                 return;\r
 112         }\r
 113         syndrome = (u8*)mxGetPr(prhs[MATLAB_in_message]);\r
 114         syndromelength = (int)mxGetM(prhs[MATLAB_in_message]);\r
 115 \r
 116         if(!mxIsDouble(prhs[MATLAB_in_pricevec])) {\r
 117                 if(!mxIsUint8(prhs[MATLAB_in_pricevec])) {\r
 118                         mexErrMsgTxt("The price vector must be of type double or uint8.");\r
 119                         return;\r
 120                 } else {\r
 121                         usefloat = false;\r
 122                 }\r
 123         } else {\r
 124                 usefloat = true;\r
 125         }\r
 126         pricevectorv = (void*)mxGetPr(prhs[MATLAB_in_pricevec]);\r
 127 \r
 128         if(nrhs > 3)\r
 129                 matrixheight = (int)mxGetScalar(prhs[MATLAB_in_constrLength]);\r
 130         else\r
 131                 matrixheight = 10;\r
 132 \r
 133         // end of MATLAB interface\r
 134 \r
 135         if(matrixheight > 31) {\r
 136                 mexErrMsgTxt("Submatrix height must not exceed 31.");\r
 137                 return;\r
 138         }\r
 139 \r
 140         height = 1 << matrixheight;\r
 141         colmask = height - 1;\r
 142         height = (height + 31) & (~31);\r
 143 \r
 144         parts = height >> 5;\r
 145         \r
 146         if(nlhs > 1) {\r
 147                 path = (u32*)malloc(vectorlength * parts * sizeof(u32));\r
 148                 if(path == NULL) {\r
 149                         char errmsg[100];\r
 150                         sprintf(errmsg, "Not enough memory (%u byte array could not be allocated).\n", vectorlength * parts * sizeof(u32));\r
 151                         mexErrMsgTxt(errmsg);\r
 152                         return;\r
 153                 }\r
 154                 pathindex = 0;\r
 155         }\r
 156 \r
 157         {\r
 158                 int shorter, longer, worm;\r
 159                 double invalpha;\r
 160 \r
 161                 matrices = (int *)malloc(syndromelength * sizeof(int));\r
 162                 widths = (int *)malloc(syndromelength * sizeof(int));\r
 163 \r
 164                 invalpha = (double)vectorlength / syndromelength;\r
 165                 if(invalpha < 1) {\r
 166                         free(matrices);\r
 167                         free(widths);\r
 168                         if(nlhs > 1)\r
 169                                 free(path);\r
 170                         mexErrMsgTxt("The message cannot be longer than the cover object.\n");\r
 171                         return;\r
 172                 }\r
 173                 if(invalpha < 2) {\r
 174                         mexWarnMsgTxt("The relative payload is greater than 1/2. This may result in poor embedding efficiency.\n");\r
 175                 }\r
 176                 shorter = (int)floor(invalpha);\r
 177                 longer = (int)ceil(invalpha);\r
 178                 if((columns[0] = getMatrix(shorter, matrixheight)) == NULL) {\r
 179                         free(matrices);\r
 180                         free(widths);\r
 181                         if(nlhs > 1)\r
 182                                 free(path);\r
 183                         return;\r
 184                 }\r
 185                 if((columns[1] = getMatrix(longer, matrixheight)) == NULL) {\r
 186                         free(columns[0]);\r
 187                         free(matrices);\r
 188                         free(widths);\r
 189                         if(nlhs > 1)\r
 190                                 free(path);\r
 191                         return;\r
 192                 }\r
 193                 worm = 0;\r
 194                 for(i = 0; i < syndromelength; i++) {\r
 195                         if(worm + longer <= (i + 1) * invalpha + 0.5) {\r
 196                                 matrices[i] = 1;\r
 197                                 widths[i] = longer;\r
 198                                 worm += longer;\r
 199                         } else {\r
 200                                 matrices[i] = 0;\r
 201                                 widths[i] = shorter;\r
 202                                 worm += shorter;\r
 203                         }\r
 204                 }\r
 205         }\r
 206 \r
 207         if(usefloat) {\r
 208         /*\r
 209                 SSE FLOAT VERSION\r
 210         */\r
 211                 int pathindex8 = 0;\r
 212                 int shift[2] = {0, 4};\r
 213                 u8 mask[2] = {0xf0, 0x0f};\r
 214                 float *prices;\r
 215                 u8 *path8 = (u8*)path;\r
 216                 double *pricevector = (double*)pricevectorv;\r
 217                 double total = 0;\r
 218                 double inf = std::numeric_limits<double>::infinity();\r
 219 \r
 220                 sseheight = height >> 2;\r
 221                 ssedone = (u8*)malloc(sseheight * sizeof(u8));\r
 222                 prices = (float*)aligned_malloc(height * sizeof(float), 16);\r
 223 \r
 224                 {\r
 225                         __m128 fillval = _mm_set1_ps(inf);\r
 226                         for(i = 0; i < height; i+= 4) {\r
 227                                 _mm_store_ps(&prices[i], fillval);\r
 228                                 ssedone[i >> 2] = 0;\r
 229                         }\r
 230                 }\r
 231 \r
 232                 prices[0] = 0.0f;\r
 233 \r
 234                 for(index = 0, index2 = 0; index2 < syndromelength; index2++) {\r
 235                         register __m128 c1, c2;\r
 236 \r
 237                         for(k = 0; k < widths[index2]; k++, index++) {\r
 238                                 column = columns[matrices[index2]][k] & colmask;\r
 239 \r
 240                                 if(vector[index] == 0) {\r
 241                                         c1 = _mm_setzero_ps();\r
 242                                         c2 = _mm_set1_ps((float)pricevector[index]);\r
 243                                 } else {\r
 244                                         c1 = _mm_set1_ps((float)pricevector[index]);\r
 245                                         c2 = _mm_setzero_ps();\r
 246                                 }\r
 247 \r
 248                                 total += pricevector[index];\r
 249 \r
 250                                 for(m = 0; m < sseheight; m++) {\r
 251                                         if(!ssedone[m]) {\r
 252                                                 register __m128 v1, v2, v3, v4;\r
 253                                                 altm = (m ^ (column >> 2));\r
 254                                                 v1 = _mm_load_ps(&prices[m << 2]);\r
 255                                                 v2 = _mm_load_ps(&prices[altm << 2]);\r
 256                                                 v3 = v1;\r
 257                                                 v4 = v2;\r
 258                                                 ssedone[m] = 1;\r
 259                                                 ssedone[altm] = 1;\r
 260                                                 switch(column & 3) {\r
 261                                                 case 0:\r
 262                                                         break;\r
 263                                                 case 1:\r
 264                                                         v2 = _mm_shuffle_ps(v2, v2, 0xb1);\r
 265                                                         v3 = _mm_shuffle_ps(v3, v3, 0xb1);\r
 266                                                         break;\r
 267                                                 case 2:\r
 268                                                         v2 = _mm_shuffle_ps(v2, v2, 0x4e);\r
 269                                                         v3 = _mm_shuffle_ps(v3, v3, 0x4e);\r
 270                                                         break;\r
 271                                                 case 3:\r
 272                                                         v2 = _mm_shuffle_ps(v2, v2, 0x1b);\r
 273                                                         v3 = _mm_shuffle_ps(v3, v3, 0x1b);\r
 274                                                         break;\r
 275                                                 }\r
 276                                                 v1 = _mm_add_ps(v1, c1);\r
 277                                                 v2 = _mm_add_ps(v2, c2);\r
 278                                                 v3 = _mm_add_ps(v3, c2);\r
 279                                                 v4 = _mm_add_ps(v4, c1);\r
 280 \r
 281                                                 v1 = _mm_min_ps(v1, v2);\r
 282                                                 v4 = _mm_min_ps(v3, v4);\r
 283                                                                                         \r
 284                                                 _mm_store_ps(&prices[m << 2], v1);\r
 285                                                 _mm_store_ps(&prices[altm << 2], v4);\r
 286 \r
 287                                                 if(nlhs > 1) {\r
 288                                                         v2 = _mm_cmpeq_ps(v1, v2);\r
 289                                                         v3 = _mm_cmpeq_ps(v3, v4);\r
 290                                                         path8[pathindex8 + (m >> 1)] = (path8[pathindex8 + (m >> 1)] & mask[m & 1]) | (_mm_movemask_ps(v2) << shift[m & 1]);\r
 291                                                         path8[pathindex8 + (altm >> 1)] = (path8[pathindex8 + (altm >> 1)] & mask[altm & 1]) | (_mm_movemask_ps(v3) << shift[altm & 1]);\r
 292                                                 }\r
 293                                         }\r
 294                                 }\r
 295 \r
 296                                 for(i = 0; i < sseheight; i++) {\r
 297                                         ssedone[i] = 0;\r
 298                                 }\r
 299 \r
 300                                 pathindex += parts;\r
 301                                 pathindex8 += parts << 2;\r
 302                         }\r
 303 \r
 304                         if(syndrome[index2] == 0) {\r
 305                                 for(i = 0, l = 0; i < sseheight; i+= 2, l += 4) {\r
 306                                         _mm_store_ps(&prices[l], _mm_shuffle_ps(_mm_load_ps(&prices[i << 2]), _mm_load_ps(&prices[(i + 1) << 2]), 0x88));\r
 307                                 }\r
 308                         } else {\r
 309                                 for(i = 0, l = 0; i < sseheight; i+= 2, l += 4) {\r
 310                                         _mm_store_ps(&prices[l], _mm_shuffle_ps(_mm_load_ps(&prices[i << 2]), _mm_load_ps(&prices[(i + 1) << 2]), 0xdd));\r
 311                                 }\r
 312                         }\r
 313 \r
 314                         if(syndromelength - index2 <= matrixheight)\r
 315                                 colmask >>= 1;\r
 316 \r
 317                         {\r
 318                                 register __m128 fillval = _mm_set1_ps(inf);\r
 319                                 for(l >>= 2; l < sseheight; l++) {\r
 320                                         _mm_store_ps(&prices[l << 2], fillval);\r
 321                                 }\r
 322                         }\r
 323                 }\r
 324 \r
 325                 totalprice = prices[0];\r
 326 \r
 327                 aligned_free(prices);\r
 328                 free(ssedone);\r
 329 \r
 330                 if(totalprice >= total) {\r
 331                         free(matrices);\r
 332                         free(widths);\r
 333                         free(columns[0]);\r
 334                         free(columns[1]);\r
 335                         if(nlhs > 1)\r
 336                                 free(path);\r
 337                         mexErrMsgTxt("The syndrome is not in the range of the syndrome matrix.\n");\r
 338                         return;\r
 339                 }\r
 340         } else {\r
 341         /*\r
 342                 SSE UINT8 VERSION\r
 343         */\r
 344                 int pathindex16 = 0, subprice = 0;\r
 345                 u8 maxc = 0, minc = 0;\r
 346                 u8 *prices, *pricevector = (u8*)pricevectorv;\r
 347                 u16 *path16 = (u16 *)path;\r
 348                 __m128i *prices16B;\r
 349 \r
 350                 sseheight = height >> 4;\r
 351                 ssedone = (u8*)malloc(sseheight * sizeof(u8));\r
 352                 prices = (u8*)aligned_malloc(height * sizeof(u8), 16);\r
 353                 prices16B = (__m128i*)prices;\r
 354 \r
 355                 {\r
 356                         __m128i napln = _mm_set1_epi32(0xffffffff);\r
 357                         for(i = 0; i < sseheight; i++) {\r
 358                                 _mm_store_si128(&prices16B[i], napln);\r
 359                                 ssedone[i] = 0;\r
 360                         }\r
 361                 }\r
 362 \r
 363                 prices[0] = 0;\r
 364 \r
 365                 for(index = 0, index2 = 0; index2 < syndromelength; index2++) {\r
 366                         register __m128i c1, c2, maxp, minp;\r
 367 \r
 368                         if((u32)maxc + pricevector[index] >= 254) {\r
 369                                 aligned_free(path);\r
 370                                 free(ssedone);\r
 371                                 free(matrices);\r
 372                                 free(widths);\r
 373                                 free(columns[0]);\r
 374                                 free(columns[1]);\r
 375                                 if(nlhs > 1)\r
 376                                         free(path);\r
 377                                 mexErrMsgTxt("Price vector limit exceeded.");\r
 378                                 return;\r
 379                         }\r
 380 \r
 381                         for(k = 0; k < widths[index2]; k++, index++) {\r
 382                                 column = columns[matrices[index2]][k] & colmask;\r
 383 \r
 384                                 if(vector[index] == 0) {\r
 385                                         c1 = _mm_setzero_si128();\r
 386                                         c2 = _mm_set1_epi8(pricevector[index]);\r
 387                                 } else {\r
 388                                         c1 = _mm_set1_epi8(pricevector[index]);\r
 389                                         c2 = _mm_setzero_si128();\r
 390                                 }\r
 391 \r
 392                                 minp = _mm_set1_epi8(-1);\r
 393                                 maxp = _mm_setzero_si128();\r
 394 \r
 395                                 for(m = 0; m < sseheight; m++) {\r
 396                                         if(!ssedone[m]) {\r
 397                                                 register __m128i v1, v2, v3, v4;\r
 398                                                 altm = (m ^ (column >> 4));\r
 399                                                 v1 = _mm_load_si128(&prices16B[m]);\r
 400                                                 v2 = _mm_load_si128(&prices16B[altm]);\r
 401                                                 v3 = v1;\r
 402                                                 v4 = v2;\r
 403                                                 ssedone[m] = 1;\r
 404                                                 ssedone[altm] = 1;\r
 405                                                 if(column & 8) {\r
 406                                                         v2 = _mm_shuffle_epi32(v2, 0x4e);\r
 407                                                         v3 = _mm_shuffle_epi32(v3, 0x4e);\r
 408                                                 }\r
 409                                                 if(column & 4) {\r
 410                                                         v2 = _mm_shuffle_epi32(v2, 0xb1);\r
 411                                                         v3 = _mm_shuffle_epi32(v3, 0xb1);\r
 412                                                 }\r
 413                                                 if(column & 2) {\r
 414                                                         v2 = _mm_shufflehi_epi16(v2, 0xb1);\r
 415                                                         v3 = _mm_shufflehi_epi16(v3, 0xb1);\r
 416                                                         v2 = _mm_shufflelo_epi16(v2, 0xb1);\r
 417                                                         v3 = _mm_shufflelo_epi16(v3, 0xb1);\r
 418                                                 }\r
 419                                                 if(column & 1) {\r
 420                                                         v2 = _mm_or_si128(_mm_srli_epi16(v2, 8), _mm_slli_epi16(v2, 8));\r
 421                                                         v3 = _mm_or_si128(_mm_srli_epi16(v3, 8), _mm_slli_epi16(v3, 8));\r
 422                                                 }\r
 423                                                 v1 = _mm_adds_epu8(v1, c1);\r
 424                                                 v2 = _mm_adds_epu8(v2, c2);\r
 425                                                 v3 = _mm_adds_epu8(v3, c2);\r
 426                                                 v4 = _mm_adds_epu8(v4, c1);\r
 427 \r
 428                                                 v1 = _mm_min_epu8(v1, v2);\r
 429                                                 v4 = _mm_min_epu8(v3, v4);\r
 430                                                                                         \r
 431                                                 _mm_store_si128(&prices16B[m], v1);\r
 432                                                 _mm_store_si128(&prices16B[altm], v4);\r
 433 \r
 434                                                 minp = _mm_min_epu8(minp, _mm_min_epu8(v1, v4));\r
 435                                                 maxp = _mm_max_epu8(maxp, maxLessThan255(v1, v4));\r
 436                                                 \r
 437                                                 if(nlhs > 1) {\r
 438                                                         v2 = _mm_cmpeq_epi8(v1, v2);\r
 439                                                         v3 = _mm_cmpeq_epi8(v3, v4);\r
 440                                                         path16[pathindex16 + m] = (u16)_mm_movemask_epi8(v2);\r
 441                                                         path16[pathindex16 + altm] = (u16)_mm_movemask_epi8(v3);\r
 442                                                 }\r
 443                                         }\r
 444                                 }\r
 445 \r
 446                                 maxc = max16B(maxp);\r
 447                                 minc = min16B(minp);\r
 448 \r
 449                                 maxc -= minc;\r
 450                                 subprice += minc;\r
 451                                 {\r
 452                                         register __m128i mask = _mm_set1_epi32(0xffffffff);\r
 453                                         register __m128i m = _mm_set1_epi8(minc);\r
 454                                         for(i = 0; i < sseheight; i++) {\r
 455                                                 register __m128i res;\r
 456                                                 register __m128i pr = prices16B[i];\r
 457                                                 res = _mm_andnot_si128(_mm_cmpeq_epi8(pr, mask), m);\r
 458                                                 prices16B[i] = _mm_sub_epi8(pr, res);\r
 459                                                 ssedone[i] = 0;\r
 460                                         }\r
 461                                 }\r
 462 \r
 463                                 pathindex += parts;\r
 464                                 pathindex16 += parts << 1;\r
 465                         }\r
 466 \r
 467                         {\r
 468                                 register __m128i mask = _mm_set1_epi32(0x00ff00ff);\r
 469 \r
 470                                 if(minc == 255) {\r
 471                                         aligned_free(path);\r
 472                                         free(ssedone);\r
 473                                         free(matrices);\r
 474                                         free(widths);\r
 475                                         free(columns[0]);\r
 476                                         free(columns[1]);\r
 477                                         if(nlhs > 1)\r
 478                                                 free(path);\r
 479                                         mexErrMsgTxt("The syndrome is not in the syndrome matrix range.\n");\r
 480                                         return;\r
 481                                 }\r
 482 \r
 483                                 if(syndrome[index2] == 0) {\r
 484                                         for(i = 0, l = 0; i < sseheight; i += 2, l++) {\r
 485                                                 _mm_store_si128(&prices16B[l], _mm_packus_epi16(_mm_and_si128(_mm_load_si128(&prices16B[i]), mask), \r
 486                                                         _mm_and_si128(_mm_load_si128(&prices16B[i + 1]), mask)));\r
 487                                         }\r
 488                                 } else {\r
 489                                         for(i = 0, l = 0; i < sseheight; i += 2, l++) {\r
 490                                                 _mm_store_si128(&prices16B[l], _mm_packus_epi16(_mm_and_si128(_mm_srli_si128(_mm_load_si128(&prices16B[i]), 1), mask), \r
 491                                                         _mm_and_si128(_mm_srli_si128(_mm_load_si128(&prices16B[i + 1]), 1), mask)));\r
 492                                         }\r
 493                                 }\r
 494 \r
 495                                 if(syndromelength - index2 <= matrixheight)\r
 496                                         colmask >>= 1;\r
 497 \r
 498                                 register __m128i fillval = _mm_set1_epi32(0xffffffff);\r
 499                                 for(; l < sseheight; l++)\r
 500                                         _mm_store_si128(&prices16B[l], fillval);\r
 501                         }\r
 502                 }\r
 503                         \r
 504                 totalprice = subprice + prices[0];\r
 505                 \r
 506                 aligned_free(prices);\r
 507                 free(ssedone);\r
 508         }\r
 509         \r
 510         if(nlhs > 1) {\r
 511                 int size[2];\r
 512                 u8 *closestData;\r
 513                 mxArray *closest;\r
 514 \r
 515                 size[0] = vectorlength;\r
 516                 size[1] = 1;\r
 517                 \r
 518                 closest = mxCreateNumericArray(2, size, mxUINT8_CLASS, mxREAL);\r
 519                 closestData = (u8*)mxGetPr(closest);\r
 520 \r
 521                 pathindex -= parts;\r
 522                 index--;\r
 523                 index2--;\r
 524                 state = 0;\r
 525 \r
 526                 int h = syndromelength;\r
 527                 state = 0;\r
 528                 colmask = 0;\r
 529                 for(; index2 >= 0; index2--) {\r
 530                         for(k = widths[index2] - 1; k >= 0; k--, index--) {\r
 531                                 if(k == widths[index2] - 1) {\r
 532                                         state = (state << 1) | syndrome[index2];\r
 533                                         if(syndromelength - index2 <= matrixheight)\r
 534                                                 colmask = (colmask << 1) | 1;\r
 535                                 }\r
 536 \r
 537                                 if(path[pathindex + (state >> 5)] & (1 << (state & 31))) {\r
 538                                         closestData[index] = 1;\r
 539                                         state = state ^ (columns[matrices[index2]][k] & colmask);\r
 540                                 } else {\r
 541                                         closestData[index] = 0;\r
 542                                 }\r
 543 \r
 544                                 pathindex -= parts;\r
 545                         }\r
 546                 }\r
 547                 plhs[MATLAB_out_stego] = closest;\r
 548                 free(path);\r
 549         }\r
 550 \r
 551         free(matrices);\r
 552         free(widths);\r
 553         free(columns[0]);\r
 554         free(columns[1]);\r
 555 \r
 556         if(nlhs > 0) {\r
 557                 int size[2];\r
 558                 mxArray *arr;\r
 559 \r
 560                 size[0] = size[1] = 1;\r
 561                 arr = mxCreateNumericArray(2, size, mxDOUBLE_CLASS, mxREAL);\r
 562                 (mxGetPr(arr))[0] = totalprice;\r
 563                 plhs[MATLAB_out_price] = arr;\r
 564         }\r
 565 \r
 566         return;\r
 567 }\r