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