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