code/canny_p3-p7-stc.py

   1 #-*- coding:utf-8 -*-
   2 import Image as im
   3 import numpy
   4 from sys import argv
   5 from pylab import *
   6 import string
   7 import cv
   8 import os
   9 from random import *
  10 from math import *
  11
  12 infinity = 1000000000
  13
  14
  15 #paths for cover and stego
  16 path_stego = '/home/couturie/ajeter/stego/'
  17 path_cover = '/home/couturie/ajeter/cover/'
  18
  19 # forward part
  20
  21
  22 def dec(ch,n):
  23     l = len(ch)
  24     acc = 0
  25     for i in xrange(l):
  26         if ch[i]==1:
  27             acc = acc + 2**(n-i-1)
  28     return acc
  29
  30
  31 def bin(elem,n):
  32     """Convertit un nombre en binaire"""
  33     q = -1
  34     res = [0 for i in xrange(n)]
  35     i = 1
  36     while q != 0:
  37         q = elem // 2
  38         r = elem % 2
  39         res[n-i] =  r
  40         elem = q
  41         i+=1
  42     return res
  43
  44
  45
  46 def xorb(a,b):
  47     return 1 if a != b else 0
  48
  49 def xor(e1,e2,h):
  50     e1b,e2b  = bin(e1,h),bin(e2,h)
  51     d = dec([xorb(e1b[j],e2b[j]) for j in xrange(h)],h)
  52     return d
  53
  54 def lit(d,(indx,indy)):
  55     if (indx,indy) in d :
  56         return d[(indx,indy)]
  57     else :
  58         return 0
  59
  60
  61
  62
  63 def forward(H_hat,x,message,lnm,rho):
  64     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
  65     path = dict()
  66     nbblock = lnm
  67     wght = [infinity for _ in xrange(int(2**h))]
  68     wght[0]=0
  69     newwght = [0 for _ in xrange(int(2**h))]
  70 #    rho = 1
  71 #    rho= [1 for _ in xrange(len(x))]
  72     indx,indm = 0,0
  73     i=0
  74     while i < nbblock: # pour chaque bit du message
  75         for j in xrange(w):   # pour chaque colonne de H_hat
  76             #print indx, "en entrant",wght
  77             k = 0
  78             while k < int(2**h): # pour chaque ligne de H
  79                 w0 = wght[k] + x[indx]*rho[indx]
  80                 w1 = wght[xor(k,H_hat[j],h)] + (1-x[indx])*rho[indx]
  81                 if w1 < w0 :
  82                     path[(indx,k)] = 1
  83                 else :
  84                     if (indx,k) in path:
  85                         del path[(indx,k)]
  86                 newwght[k] = min(w0,w1)
  87                 k +=1
  88             indx +=1
  89             wght = [t for t in newwght]
  90             #print " apres calcul",wght
  91
  92         for j in xrange(int(2**(h-1))):   # pour chaque colonne de H
  93             wght[j] = wght[2*j + message[indm]]
  94         wght = wght[:int(pow(2,h-1))] + [infinity for _ in xrange(int(pow(2,h)-pow(2,h-1)))]
  95         indm +=1
  96         i +=1
  97     start = np.argmin(wght)
  98     return (start,path)
  99
 100
 101 def backward(start,H_hat,x,message,lnm,path):
 102     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 103     indx,indm = len(x)-1,lnm-1
 104     state = 2*start + message[indm]
 105     indm -=1
 106     # l'initialisation de state n'est pas optimale...
 107     nbblock = lnm
 108     y=np.zeros(len(x))
 109     i=0
 110     while i < nbblock:
 111         l = range(w)
 112         l.reverse()
 113         for j in l:   # pour chaque colonne de H_hat
 114             y[indx] = lit(path,(indx,state))
 115             state = xor(state,y[indx]*H_hat[j],h)
 116             indx -=1
 117         state = 2*state + message[indm]
 118         indm -=1
 119         i +=1
 120     return [int(t) for t in y]
 121
 122
 123
 124
 125
 126
 127 def trouve_H_hat(n,m,h):
 128     assert h ==7
 129     alpha = float(n)/m
 130     assert alpha >= 1
 131     index = min(int(alpha),9)
 132     mat = {
 133         1 : [127],
 134         2 : [71,109],
 135         3 : [95, 101, 121],
 136         4 : [81, 95, 107, 121],
 137         5 : [75, 95, 97, 105, 117],
 138         6 : [73, 83, 95, 103, 109, 123],
 139         7 : [69, 77, 93, 107, 111, 115, 121],
 140         8 : [69, 79, 81, 89, 93, 99, 107, 119],
 141         9 : [69, 79, 81, 89, 93, 99, 107, 119, 125]
 142         }[index]
 143     return(mat, index*m)
 144
 145
 146 def stc(x,rho,message):
 147     lnm = len(message)
 148     (mat,taille_suff) = trouve_H_hat(len(x),len(message),7)
 149     x_b = x[:taille_suff]
 150     (start,path) = forward(mat,x_b,message,lnm,rho)
 151     return (x_b,backward(start,mat,x_b,message,lnm,path),mat)
 152
 153
 154
 155
 156
 157 def nbdif(x,y):
 158     r,it = 0,0
 159     l = len(y)
 160     while it < l :
 161         if x[it] != y[it] :
 162             r +=1
 163         it += 1
 164     return float(r)/l
 165
 166
 167
 168
 169
 170
 171 def prod(H_hat,lnm,y):
 172     (h,w) = int(log(max(H_hat),2))+1, len(H_hat)
 173     i=0
 174     H =[]
 175     V=[0 for _ in range(len(y))]
 176     sol=[]
 177     while i < lnm: # pour chaque ligne
 178         V=[0 for _ in range(len(y))]
 179         k = max([(i-h+1)*w,0])
 180         dec = max([i-h+1,0])
 181         for j in xrange(min([i+1,h])): #nbre de blocks presents sur la ligne i
 182             for l in xrange(w): # pour chaque collone de H_hat
 183                 V[k] = bin(H_hat[l],h)[h-i-1+j+dec]
 184                 k+=1
 185
 186         sol.append(np.dot(np.array(V),np.array(y)))
 187         i+=1
 188         #H += [V]
 189     #H = np.array(H)
 190     #y = np.array(y)
 191     #print "dot",np.dot(H,y),H.shape
 192     #print "sol",sol
 193     return sol#list(np.dot(H,y))
 194
 195 def equiv(x,y):
 196     lx = len(x)
 197     assert lx == len(y)
 198     i=0
 199     while i < lx :
 200         if x[i] % 2 != y[i]%2 :
 201             return False
 202         i += 1
 203     return True
 204
 205
 206
 207
 208
 209
 210 def conversion(nombre, base, epsilon = 0.00001 ):
 211     ''' Soit nombre écrit en base 10, le retourne en base base'''
 212     if not 2 <= base <= 36:
 213         raise ValueError, "La base doit être entre 2 et 36"
 214     if not base == 2 and '.' in str(nombre):
 215         raise ValueError, "La partie décimale n'est pas gérée pour les bases\
 216                            différentes de 2."
 217     # IMPROVE : Convertir aussi la partie décimale, quand la base n'est pas égale
 218     # à 2.
 219     abc = string.digits + string.letters
 220     result = ''
 221     if nombre < 0:
 222         nombre = -nombre
 223         sign = '-'
 224     else:
 225         sign = ''
 226     if '.' in str(nombre):
 227         entier,decimal = int(str(nombre).split('.')[0]),\
 228                          float('.'+str(nombre).split('.')[1])
 229     else:
 230         entier,decimal = int(str(nombre)),0
 231     while entier !=0 :
 232         entier, rdigit = divmod( entier, base )
 233         result = abc[rdigit] + result
 234     flotante, decimalBin = 1./float(base),''
 235     while flotante > epsilon :
 236         if decimal >= flotante:
 237             decimalBin+='1'
 238             decimal-=flotante
 239         else :
 240             decimalBin+='0'
 241         flotante = flotante/float(base)
 242     if '1' in decimalBin :
 243         reste = '.'+decimalBin
 244         while reste[-1]=='0':
 245             reste = reste[:-1]
 246     else :
 247         reste = ''
 248     return sign + result + reste
 249
 250
 251 def getBit(X,pos):
 252     '''Récupère le bit en position pos de X.
 253     Par exemple, getBit(8,1) = 0, puisque le bit le plus à droite de 8 = 1000 est 0.
 254     On fera attention à ce que :
 255         - on compte à partir du point,
 256         - l'élément juste à gauche du point est en position 1,
 257         - celui juste à droite est en position -1.'''
 258     assert pos != 0
 259     entier = conversion(X,2)
 260     if '.' in entier:
 261         entier, decimal = entier.split('.')
 262         if decimal == '0':
 263             decimal = ''
 264     else:
 265         decimal = ''
 266     if '-' in entier:
 267         entier = entier.replace('-','')
 268     entier  = entier.zfill(abs(pos))
 269     decimal = (decimal+'0'*abs(pos))[:max(len(decimal),abs(pos))]
 270
 271     return int(entier[len(entier)-pos]) if pos >0 else int(decimal[-pos-1])
 272
 273
 274 def setBit(X,pos,y):
 275     '''Fixe le bit pos de X à la valeur y.
 276     Le fonctionnement est similaire à getBit :
 277         - on compte à partir du point,
 278         - l'élément juste à gauche du point est en position 1,
 279         - celui juste à droite est en position -1.'''
 280     assert pos != 0
 281     entier = conversion(X,2)
 282     if '.' in entier:
 283         entier, decimal = entier.split('.')
 284     else:
 285         decimal = ''
 286     entier  = list(entier.zfill(abs(pos)))
 287     decimal = list((decimal+'0'*abs(pos))[:max(len(decimal),abs(pos))])
 288     if pos>0:
 289         entier[len(entier)-pos]=str(int(y))
 290     else:
 291         decimal[-pos-1] = str(int(y))
 292     if decimal == []:
 293         return int(''.join(entier),2)
 294     else:
 295         S=0
 296         for k in range(len(decimal)):
 297             S += 1./2**(k+1)*int(decimal[k])
 298         return float(str(int(''.join(entier),2))+'.'+str(S).split('.')[1])
 299
 300
 301 def a2b(a):
 302     ai = ord(a)
 303     return ''.join('01'[(ai >> x) & 1] for x in xrange(7, -1, -1))
 304
 305
 306
 307 def a2b_list(L):
 308     LL=[]
 309     for i in L:
 310         for j in list(a2b(i)):
 311             LL.append(j)
 312     return LL
 313
 314
 315
 316 def toDecimal(x):
 317     return sum(map(lambda z: int(x[z]) and 2**(len(x) - z - 1),
 318                    range(len(x)-1, -1, -1)))
 319
 320 def conv_list_bit(L):
 321     L2=[]
 322     for j in range(len(L)/8):
 323         L2.append(chr(toDecimal("".join(L[j*8:(j+1)*8]))))
 324     return ''.join(L2)
 325
 326 def Denary2Binary(n):
 327     '''convert denary integer n to binary string bStr'''
 328     bStr = ''
 329     if n < 0:  raise ValueError, "must be a positive integer"
 330     if n == 0: return '0'
 331     while n > 0:
 332         bStr = str(n % 2) + bStr
 333         n = n >> 1
 334     return bStr
 335
 336
 337 def compute_filter_sobel(level,image):
 338     level2=level.copy()
 339     level2= array(level2.getdata()).flatten()
 340     l=0
 341     for x in level2:
 342         level2[l]=(x/2)*2
 343         l+=1
 344     level2_im=im.new('L',image.size)
 345     level2_im.putdata(level2)
 346
 347     cv_im = cv.CreateImageHeader(image.size, cv.IPL_DEPTH_8U, 1)
 348     cv.SetData(cv_im, level2_im.tostring())
 349     dst16 = cv.CreateImage(cv.GetSize(cv_im), cv.IPL_DEPTH_16S, 1)
 350
 351     laplace = cv.Sobel(cv_im, dst16,1, 1,7)
 352
 353     dst8 = cv.CreateImage (cv.GetSize(cv_im), cv.IPL_DEPTH_8U, 1)
 354     cv.ConvertScale(dst16,dst8)
 355     processed=im.fromstring("L", cv.GetSize(dst8), dst8.tostring())
 356 #    cv.ShowImage ('canny', dst8)
 357 #    cv.WaitKey()
 358
 359     return processed
 360
 361
 362
 363 def compute_list_bit_to_change(threshold,processed):
 364     List=[]
 365     nb=0
 366     l=0
 367     for i in processed:
 368         if (processed[l]>=threshold):
 369             #if nb%2==0:
 370                 List.append(l)
 371             #nb+=1
 372         l+=1
 373     return List
 374
 375
 376 def compute_filter_canny(level,image):
 377     level2=level.copy()
 378     level2= array(level2.getdata()).flatten()
 379     l=0
 380     for x in level2:
 381         level2[l]=(x/2)*2
 382         l+=1
 383     level2_im=im.new('L',image.size)
 384     level2_im.putdata(level2)
 385     level2_im=im.merge("RGB",(level2_im,level2_im,level2_im))
 386
 387     mean=numpy.mean(level2)
 388     std=numpy.std(level2)
 389
 390     cv_im = cv.CreateImageHeader(image.size, cv.IPL_DEPTH_8U, 3)
 391     cv.SetData(cv_im, level2_im.tostring())
 392
 393     yuv = cv.CreateImage(cv.GetSize(cv_im), 8, 3)
 394     gray = cv.CreateImage(cv.GetSize(cv_im), 8, 1)
 395     cv.CvtColor(cv_im, yuv, cv.CV_BGR2YCrCb)
 396     cv.Split(yuv, gray, None, None, None)
 397
 398     canny = cv.CreateImage(cv.GetSize(cv_im), 8, 1)
 399
 400     List_bit_to_change=set([])
 401     Weight=[]
 402
 403
 404
 405     cv.Canny(gray, canny, mean-1*std, mean+1*std,3)  #avant 10 255
 406     processed=im.fromstring("L", cv.GetSize(canny), canny.tostring())
 407     processed= array(processed.getdata()).flatten()
 408     List3=set(compute_list_bit_to_change(100,processed))
 409
 410     cv.Canny(gray, canny, mean-1*std, mean+1*std,5)  #avant 10 255
 411     processed=im.fromstring("L", cv.GetSize(canny), canny.tostring())
 412     processed= array(processed.getdata()).flatten()
 413     List5=set(compute_list_bit_to_change(100,processed))
 414
 415     cv.Canny(gray, canny, mean-1*std, mean+1*std,7)  #avant 10 255
 416     processed=im.fromstring("L", cv.GetSize(canny), canny.tostring())
 417     processed= array(processed.getdata()).flatten()
 418     List7=set(compute_list_bit_to_change(100,processed))
 419
 420     nb_bit_embedded=(512*512/10)+40
 421     AvailablePixel3=List3
 422     AvailablePixel5=AvailablePixel3.union(List5)
 423     AvailablePixel7=AvailablePixel5.union(List7)
 424     if len(AvailablePixel3)>nb_bit_embedded:
 425         step=1
 426         WorkingPixel=AvailablePixel3
 427     elif len(AvailablePixel5)>nb_bit_embedded:
 428         step=2
 429         WorkingPixel=AvailablePixel5
 430     elif len(AvailablePixel7)>nb_bit_embedded:
 431         step=3
 432         WorkingPixel=AvailablePixel7
 433     else:
 434         step=4
 435         WorkingPixel=range(len(level2))
 436
 437     print "avail P3",len(AvailablePixel3)
 438     print "avail P5",len(AvailablePixel5)
 439     print "avail P7",len(AvailablePixel7)
 440
 441     print "size WorkingPixel",len(WorkingPixel)
 442     Weight=[0 for _ in WorkingPixel]
 443
 444     l=0
 445     for i in WorkingPixel:
 446         if step>=1 and i in List3:
 447             Weight[l]=1
 448         if step>=2 and i in List5 and Weight[l]==0:
 449             Weight[l]=10
 450         if step>=3 and i in List7 and Weight[l]==0:
 451             Weight[l]=100
 452         if step>=4 and Weight[l]==0:
 453             Weight[l]=1000
 454         l+=1
 455
 456
 457
 458
 459     List_bit_to_change=WorkingPixel
 460
 461
 462
 463
 464
 465     return [List_bit_to_change,Weight]
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477 def mystego(filein, fileout):
 478     dd = im.open(filein)
 479     dd = dd.convert('RGB')
 480     red, green, blue = dd.split()
 481     level=red.copy()
 482
 483     [List_bit_to_change,Weight]=compute_filter_canny(level,dd)
 484     level= array(level.getdata()).flatten()
 485
 486
 487     bit_to_read=1
 488
 489
 490
 491
 492
 493
 494 #parameters for BBS
 495     M=18532395500947174450709383384936679868383424444311405679463280782405796233163977*39688644836832882526173831577536117815818454437810437210221644553381995813014959
 496     X=18532395500947174450709383384936679868383424444311
 497
 498
 499
 500
 501
 502
 503     l=0
 504
 505
 506
 507     message="Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete  Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete  Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete  Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete  Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete Ce que j'écris est très original... Bref, je suis un poete   Salut christophe, arrives tu à lire ce message? normalement tu dois lire cela. Bon voici un test avec un message un peu plus long. Bon j'en rajoute pour voir. Ce que j'écris est très original... Bref, je suis un poete voila c'est la fin blablabla:-)"
 508     message=message[0:len(message)/1]
 509     leng_msg=len(message)
 510     message=message+((leng_msg+7)/8*8-leng_msg)*" "
 511     leng_msg=len(message)
 512
 513     leng='%08d'%len(message)
 514
 515
 516
 517
 518     len_leng=len(leng)
 519     leng_error=int(len_leng)
 520     leng_cor=leng
 521     List_pack=a2b_list(leng_cor)
 522
 523     List_random=[]
 524     while len(List_random)<len(List_bit_to_change):
 525         X=(X*X)%M
 526         List_random.extend(Denary2Binary(X))
 527
 528     size=0
 529     for i in range(leng_msg/8):
 530         m=message[i*8:(i+1)*8]
 531         m_cor=m
 532         m_bin=a2b_list(m_cor)
 533         size=size+len(m_bin)
 534         List_pack.extend(m_bin)
 535
 536
 537
 538
 539
 540
 541
 542
 543     Support=[getBit(level[l],bit_to_read) for l in List_bit_to_change]
 544
 545
 546     Message=[(int(List_pack[l])^int(List_random[l])) for l in xrange(len(List_pack))]
 547
 548     print "support",len(List_bit_to_change)
 549     print "message",len(Message)
 550     print "weight",len(Weight)
 551
 552     (x_b,Stc_message,H_hat) = stc(Support,Weight,Message)
 553
 554     print "pourcentage modif",nbdif(x_b,Stc_message)
 555     print "taille Stc_message",len(Stc_message)
 556
 557
 558
 559     l=0
 560     size_mesg=0
 561     val_mod=0
 562
 563
 564
 565     for l in List_bit_to_change:
 566         if(size_mesg<len(Stc_message)):
 567             b=getBit(level[l],bit_to_read)
 568             if b!=Stc_message[size_mesg]:
 569                 val_mod+=1
 570             level[l]=float64(setBit(level[l],bit_to_read,Stc_message[size_mesg]))
 571             size_mesg+=1
 572
 573     print 'size mesg',size_mesg
 574     print 'val mod',val_mod
 575     print 'len message',len(Message),len(List_pack)
 576
 577
 578
 579     zz3=im.new('L',dd.size)
 580     zz3.putdata(level)
 581
 582     zz3.save(fileout)
 583
 584
 585
 586
 587 listing = os.listdir(path_cover)
 588
 589 print listing
 590
 591 list_nb_bit=[]
 592 l=0
 593
 594
 595
 596 for infile in listing:
 597     print "current file is: " + infile, path_stego+infile
 598     mystego(path_cover+infile,path_stego+infile)
 599     l+=1
 600