RS（Regular Singular）隐写分析及实现

上篇博客介绍了LSB（最低有效位）算法。
这篇针对LSB算法介绍RS隐写分析的实现。

先介绍几个概念：
翻转函数为：

平滑度函数为：
掩码算子M为:

掩码算子可以取随机0，1但是一定要各占50%概率。

RS分析分为几个步骤：
1.对像素分组
以8位灰度图像为例，图像的像素取值范围为集合𝑃
P=｛0，……，255｝

分组可以看情况取，不一定非要是4。
2.对分组后的数据计算相关性（用平滑度函数），分别进行F1与F-1翻转后再计算相关性，进行对比得出$R_M,S_M,R_{-M},S_{-M}$的数值。

3.对图像的每个像素最低位取反，同样进行第二步操作。得出另外一组$R_M,S_M,R_{-M},S_{-M}$
设隐写率为p，第三步算出的$R_M,S_M,R_{-M},S_{-M}$为隐写率$(1-\frac{p}{2})$情况下的，记为$R_M(1-\frac{p}{2})$

4.解下面方程，取其中较小的数为解$x$，则隐写率估计为$p=\frac{x}{x-0.5}$
判断是否经过LSB隐写：

如果经过LSB隐写则对图片进行F1翻转与进行F-1翻转会对图片平滑度造成不同的影响，进而导致$R_M,S_M,R_{-M},S_{-M}$参数的不对称性，可以明显看出经过隐写。
通过F1隐写相当于在最低位增加了噪声，加水印时也相当于增加了最低位的噪声，但是F1变换引入噪声对像素的影响最多为1，原图单独引入F-1变换引入噪声对像素的影响同样最多为1。
但是，引入F-1变换时，如果加了水印的话相当于部分像素先进行了F1变换后进行了F-1变换将噪声引入像素的影响可以达到2。

对于上文LSB算法提到的第一种方法比第二种更具有RS分析的隐蔽性，可以直观理解为原图与水印图都具有是图片具有像素上的平滑性，直接二值化嵌入的话会将水印图的像素平滑性引入加水印后的图片，这样的话导致通过平滑度计算来分析隐写率较为困难。（也不是说第一种好，第一种也有缺点，比如存储二值化图片数据、直接提取最后一位信息直接就暴露了）

下面给两个代码：一种算子是随机取0，1，一种是0，1间隔。一种平滑度计算是Z字型顺序计算，一种是一行一行的计算。

# -*- coding: utf-8 -*-
"""
Created on Wed May 13 13:15:42 2020

@author: MYM
"""
import matplotlib.pyplot as plt
import matplotlib.image as pmg
import numpy as np
import random  

def rgb2gray(rgb):

    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b

    return gray
def Fp(picture,flag=0):#单个元素的翻转 输入参数为：元素，翻转标签
    if flag==1:
        if picture%2==0:
            picture= picture+1
        else:
            picture= picture-1
    elif flag==-1:
        if picture%2==0:
            picture=picture-1
        else:
            picture=picture+1
    return picture
def Fp_block(block_ori,wei,flag=0):#块翻转，输入参数为：数据块，维度，翻转标签
    block=block_ori.copy()
    if wei==2:
        if flag==1:
            block[0,0]=Fp(block[0,0],1)
            block[0,1]=Fp(block[0,1],0)
            block[1,0]=Fp(block[1,0],1)
            block[1,1]=Fp(block[1,1],0)
        elif flag==-1:
            block[0,0]=Fp(block[0,0],-1)
            block[0,1]=Fp(block[0,1],0)
            block[1,0]=Fp(block[1,0],-1)
            block[1,1]=Fp(block[1,1],0)
    elif wei==8:
        if flag==1:
            for m in range(8):
                for n in range(8):
                    if (m*8+n)%2==0:
                        block[m,n]=Fp(block[m,n],1)
                    else:
                        block[m,n]=Fp(block[m,n],0)
        elif flag==-1:
             for m in range(8):
                for n in range(8):
                    if (m*8+n)%2==0:
                        block[m,n]=Fp(block[m,n],-1)
                    else:
                        block[m,n]=Fp(block[m,n],0)
    return block
def Smoothness(picture,f=0):#平滑度函数
    rs_picture=picture.reshape(-1,1)
    for i in range(1,rs_picture.shape[0]):
        f+=abs(rs_picture[i]-rs_picture[i-1])
    return f
def Dblock(picture,flag=2):#图片分块 flag代表方块的边长
    rs_picture=picture.copy()
    if flag==2:
        BlockSize = 4;
        BlockCol = 2;
        BlockRow = 2;
        x,y=rs_picture.shape[0],rs_picture.shape[1]
        bufsize=x*y
        Block_p= np.zeros((bufsize//BlockSize,BlockRow,BlockCol));
        Blocknum=0
        for i in range(0,x//BlockRow):
            for j in range(0,y//BlockCol):
                Block_p[Blocknum,0,0]=rs_picture[BlockRow*i,BlockCol*j]
                Block_p[Blocknum,0,1]=rs_picture[BlockRow*i,BlockCol*j+1]
                Block_p[Blocknum,1,0]=rs_picture[BlockRow*i+1,BlockCol*j]
                Block_p[Blocknum,1,1]=rs_picture[BlockRow*i+1,BlockCol*j+1]
                Blocknum = Blocknum + 1;
        return Block_p
    elif flag==8:
        BlockSize = 64;
        BlockCol = 8;
        BlockRow = 8;
        x,y=rs_picture.shape[0]//8,rs_picture.shape[1]//8
        Block_p=np.zeros((x*y,8,8))
        Blocknum=0
        for i in range(x):
            for j in range(y):
                for k in range(8):
                    for h in range(8):
                        Block_p[Blocknum,k,h]=rs_picture[i*8+k,j*8+h]
                Blocknum+=1
        return Block_p
def compu(img):#计算参数
    computer=img.copy()#需要计算的图
    RM=0
    SM=0
    NRM=0
    NSM=0
    ori=Dblock(computer,8)
    Fori=Dblock(computer,8)
    NFori=Dblock(computer,8)
    #翻转操作
    for m in range(ori.shape[0]):
        Fori[m,:,:]=Fp_block(Fori[m,:,:],8,1)
        NFori[m,:,:]=Fp_block(NFori[m,:,:],8,-1)
        if Smoothness(Fori[m,:,:])>Smoothness(ori[m,:,:]):
            RM=RM+1
        if Smoothness(Fori[m,:,:])<Smoothness(ori[m,:,:]):
            SM=SM+1
        if Smoothness(NFori[m,:,:])>Smoothness(ori[m,:,:]):
            NRM=NRM+1
        if Smoothness(NFori[m,:,:])<Smoothness(ori[m,:,:]):
            NSM=NSM+1
        
   
    return RM,SM,NRM,NSM
def LSB(fig,rate):#用来生成带水印的图
    s = int(512*512*rate)
    sec = [0] * (s)
    k = 0
    
    for i in range(s):
        sec[i] = random.randint(0,1)
    
#    print sec    
#    print fig
    
    for i in range(512):
        for j in range(512):
            if (k < s):
                if (sec[k]==1 and fig[i][j]%2==0 ):#偶数嵌入1
                    fig[i][j] = fig[i][j] + 1
                    k += 1
                elif (sec[k]==1 and fig[i][j]%2==1):#奇数嵌入1
                    fig[i][j] = fig[i][j] + 0
                    k += 1
                elif (sec[k]==0 and fig[i][j]%2==0):#偶数嵌入0
                    fig[i][j] = fig[i][j] + 0
                    k += 1
                elif (sec[k]==0 and fig[i][j]%2==1):#奇数嵌入0
                    fig[i][j] = fig[i][j] - 1
                    k += 1
    return fig


copy_right_rgb=pmg.imread("copyright.jpg")
origin_rgb=pmg.imread("lena.bmp")
#取第一个通道的图像
copy_right=np.array(copy_right_rgb[:,:,0])
origin=np.array(origin_rgb[:,:,0])
watermark=LSB(origin,0.5)
#计算RM,SM等
RM0,SM0,NRM0,NSM0=compu(watermark)
print("RM=",RM0,"SM=",SM0,"RM-=",NRM0,"SM-",NSM0)
#对图取反
watermark1=watermark.copy()
for i in range(watermark1.shape[0]):
    for j in range(watermark1.shape[1]):
        watermark1[i,j]=Fp(watermark1[i,j],1)
RM1,SM1,NRM1,NSM1=compu(watermark1)
print("RM=",RM1,"SM=",SM1,"RM-=",NRM1,"SM-",NSM1)
d0=RM0-SM0
d1=RM1-SM1
d2=NRM0-NSM0
d3=NRM1-NSM1
args=[2*(d1+d0),d2-d3-d1-3*d0,d0-d2]
a,b=np.roots(args)
if a>b:
    print("估计隐写率为：",b/(b-0.5))
else:
    print("估计隐写率为：",a/(a-0.5))

# -*- coding: utf-8 -*-
"""
Created on Fri May 15 00:55:54 2020

@author: MYM
"""


import matplotlib.pyplot as plt # plt 用于显示图片
import matplotlib.image as pmg # mpimg 用于读取图片
import numpy as np
import pylab as pl
import random  

#对像素块进行Z字排序
def Z(tmp):
    size = 8
    a = [0] * 64
    n = 0
    i = 0
    j = 0
    status = 0 #标记当前状态 0：右上运动 1：左下运动
    while (n<size*size):
        if (((i==0 and j%2!=0) or (j==1 and i==0)) and j!=size-1): #当前位于像素矩阵的最左边的奇数位且不是边界，向下移动,下一步向右上运动
            a[n] = tmp[i][j]
            j = j + 1
            n = n + 1
            status = 0    
        if (((i==0 and j%2!=0) or (j==1 and i==0)) and j==size-1): #当前位于像素矩阵的最左边的奇数位且是边界，向右移动,下一步向右上运动
            a[n] = tmp[i][j]
            i = i + 1
            n = n + 1
            status = 0    
        elif ((j==0 and i%2==0 and i>1) or (i==0 and j==0)): #当前位于像素矩阵的最上面的偶数位，向右移动，下一步向左下运动             
            a[n] = tmp[i][j]
            i = i + 1
            n = n + 1
            status = 1
        elif ((i==0 and j%2==0 and j>1) or (j==0 and i==0)): #当前位于像素矩阵的最左边的偶数位，向右上移动
            a[n] = tmp[i][j]
            j = j - 1
            i = i + 1
            n = n + 1
            status = 0
        elif ((j==0 and i%2!=0) or (i==1 and j==0)): #当前位于像素矩阵的最上面的奇数位，向左下移动
            a[n] = tmp[i][j]
            i = i - 1
            j = j + 1
            n = n + 1
            status = 1
        elif ((i==size-1 and j%2!=0)): #当前位于像素矩阵的最右边的奇数位，向下移动,下一步向左下运动
            a[n] = tmp[i][j]
            j = j + 1
            n = n + 1
            status = 1
        elif ((j==size-1 and i%2==0 and i>1)): #当前位于像素矩阵的最下面的偶数位，向右移动，下一步向右上移动              
            a[n] = tmp[i][j]
            i = i + 1
            n = n + 1
            status = 0
        elif ((i==size-1 and j%2==0 and j>1)): #当前位于像素矩阵的最右边的偶数位，向左下移动
            a[n] = tmp[i][j]
            i = i - 1
            j = j + 1
            n = n + 1
            status = 1
        elif ((j==size-1 and i%2!=0)): #当前位于像素矩阵的最下面的奇数位，向右上移动
            a[n] = tmp[i][j]
            i = i + 1
            j = j - 1
            n = n + 1
            status = 0
        else: #不是边界条件时，使用状态值判断移动方向
            if (status == 0): #右上运动
                a[n] = tmp[i][j]
                i = i + 1
                j = j - 1
                n = n + 1
                status = 0
            elif (status == 1): #左下运动
                a[n] = tmp[i][j]
                i = i - 1
                j = j + 1
                n = n + 1
                status = 1          
    return a
def LSB(fig,rate):
    s = int(512*512*rate)
    sec = [0] * (s)
    k = 0
    
    for i in range(s):
        sec[i] = random.randint(0,1)
    
#    print sec    
#    print fig
    
    for i in range(512):
        for j in range(512):
            if (k < s):
                if (sec[k]==1 and fig[i][j]%2==0 ):#偶数嵌入1
                    fig[i][j] = fig[i][j] + 1
                    k += 1
                elif (sec[k]==1 and fig[i][j]%2==1):#奇数嵌入1
                    fig[i][j] = fig[i][j] + 0
                    k += 1
                elif (sec[k]==0 and fig[i][j]%2==0):#偶数嵌入0
                    fig[i][j] = fig[i][j] + 0
                    k += 1
                elif (sec[k]==0 and fig[i][j]%2==1):#奇数嵌入0
                    fig[i][j] = fig[i][j] - 1
                    k += 1
    return fig
#计算像素相关性
def Calculate(a):
    res = 0
    for i in range(63):
        if (a[i+1] > a[i]):
            res = res + a[i+1] - a[i]
        else:
            res = res + a[i] - a[i+1]
    return res

#0翻转
def F0(val):  
    return val

#正翻转
def F1(val):
    if (val%2==0 and val!=1): #偶数加一
        val = val + 1
    elif (val%2==1 or val==1): #奇数减一
        val = val - 1
    return val
    
#负翻转
def F_1(val):
    if (val%2==0 and val!=1): #偶数减一
        val = val - 1
    elif (val%2==1 or val==1): #奇数加一
        val = val + 1
    return val
    
#生成随机数组
def Random(typ):
    ran = [0] * 64
    if (typ == 1):
        for i in range(64):
            ran[i] = random.randint(0,1)
    elif (typ == -1):
        for i in range(64):
            ran[i] = random.randint(-1,0)
    return ran

#RS隐写分析
def RS(tmp):
#==============================================================================
#     非负反转
#==============================================================================
    ran = Random(1)
    rev = [[0 for co in range(8)] for ro in range(8)] # 进行反转后的二维数组
    rm = 0
    sm = 0
    
    r1 = Z(tmp)
    res1 = Calculate(r1)#反转之前的像素相关性
    
    k = 0
    for i in range(8):
        for j in range(8):
            if (ran[k] == 0):#F0翻转
                rev[i][j] = F0(tmp[i][j])
            elif (ran[k] == 1):#F1翻转
                rev[i][j] = F1(tmp[i][j])
            k = k + 1
        
    r2 = Z(rev) # 将图像块进行Z字形排序
    res2 = Calculate(r2)#翻转之后的像素相关性
    
    if (res1 > res2):
        sm = sm + 1
    elif (res1 < res2):
        rm = rm + 1
        
#==============================================================================
#     非正翻转
#==============================================================================
    k = 0
    r_m = 0
    s_m = 0
    ran = Random(-1)
    rev = [[0 for co in range(8)] for ro in range(8)] # 进行反转后的二维数组
    for i in range(8):
        for j in range (8):
            if (ran[k] == 0):
                rev[i][j] = F0(tmp[i][j])
            elif (ran[k] == -1):
                rev[i][j] = F_1(tmp[i][j])
            k = k + 1
            
    r3 = Z(rev) # 将图像块进行Z字形排序
    res3 = Calculate(r3)#翻转之后的像素相关性
    
    if (res1 > res3):
        s_m = s_m + 1
    elif (res1 < res3):
        r_m = r_m + 1
    
    res = [rm,sm,r_m,s_m]
    return res
def Dblock(picture,flag=2):#图片分块 flag代表方块的边长
    rs_picture=picture.copy()
    if flag==2:
        BlockSize = 4;
        BlockCol = 2;
        BlockRow = 2;
        x,y=rs_picture.shape[0],rs_picture.shape[1]
        bufsize=x*y
        Block_p= np.zeros((bufsize//BlockSize,BlockRow,BlockCol));
        Blocknum=0
        for i in range(0,x//BlockRow):
            for j in range(0,y//BlockCol):
                Block_p[Blocknum,0,0]=rs_picture[BlockRow*i,BlockCol*j]
                Block_p[Blocknum,0,1]=rs_picture[BlockRow*i,BlockCol*j+1]
                Block_p[Blocknum,1,0]=rs_picture[BlockRow*i+1,BlockCol*j]
                Block_p[Blocknum,1,1]=rs_picture[BlockRow*i+1,BlockCol*j+1]
                Blocknum = Blocknum + 1;
        return Block_p
    elif flag==8:
        BlockSize = 64;
        BlockCol = 8;
        BlockRow = 8;
        x,y=rs_picture.shape[0]//8,rs_picture.shape[1]//8
        Block_p=np.zeros((x*y,8,8))
        Blocknum=0
        for i in range(x):
            for j in range(y):
                for k in range(8):
                    for h in range(8):
                        Block_p[Blocknum,k,h]=rs_picture[i*8+k,j*8+h]
                Blocknum+=1
        return Block_p


##-----main-----##
origin_rgb=pmg.imread("lena.bmp")
#取第一个通道的图像
origin=np.array(origin_rgb[:,:,0])


watermark=LSB(origin,0.9)

compu=Dblock(watermark,8)
result = [0] *4
#计算RM,SM等
for i in range(compu.shape[0]):
    res=RS(compu[i,:,:])
    for n in range(4):
                result[n] = result[n] + res[n]
print ('rm=',result[0],'  sm=',result[1],'  r_m=',result[2],'  s_m=',result[3])

#对加水印后的图像最低位取反
watermark1=watermark.copy()
for i in range(watermark1.shape[0]):
    for j in range(watermark1.shape[1]):
        watermark1[i,j]=F1(watermark1[i,j])


compu1=Dblock(watermark1,8)
result1 = [0] *4
#计算RM,SM等
for i in range(compu1.shape[0]):
    res=RS(compu1[i,:,:])
    for n in range(4):
                result1[n] = result1[n] + res[n]
print ('rm=',result1[0],'  sm=',result1[1],'  r_m=',result1[2],'  s_m=',result1[3])
#计算隐写率
d0=result[0]-result[1]
d1=result1[0]-result1[1]
d2=result[2]-result[3]
d3=result1[2]-result1[3]
args=[2*(d1+d0),d2-d3-d1-3*d0,d0-d2]
a,b=np.roots(args)
if a>b:
    print("估计隐写率为：",b/(b-0.5))
else:
    print("估计隐写率为：",a/(a-0.5))