AVERAGE函数计算平均数
函数功能

AVERAGE函数用于计算所有参数的算术平均值。

函数语法

AVERAGE(number1,number2…)

参数说明

Number1,number2,…:表示要计算平均值的1~30个参数。
实例1 忽略0值求平均分
    *案例表述*
    当需要求平均值的单元格区域中包含0值时,它们也将参与求平均值的运算。如果想排除运算区域中的0值,可以按如下方法设置公式.
    *案例实现*
    1、选中D13单元格,在公式编辑栏中输入公式:AVERAGE(IF(D2:D12<>0,D2:D12)   
    2、同时按"Ctrl+Shift+Enter"组合键,即可忽略“0”值求平均值,如图所示.



实例2     统计各班级平均分
*案例表述*
     当前表格中统计了各个班级中各学生的分数,可以利用AVERAGE函数统计出各个班级的平均分数。
*案例实现*
    1、在工作表中输入数据并建立好求解标识。
    2、选中D15单元格,在公式编辑栏中输入公式:
		=AVERAGE(IF($C$2:$C$12=$C15,D$2:D$12))
    3、同时按下”Ctrl+Shift+Enter”组合键,即可计算出“1”班的平均分数。
    4、选中D15单元格,向下复制公式到D17单元格,即可求出其他班级平均分,如图所示。



注意相对引用与绝对引用，不懂相对引用与绝对引用的百度下

哎，看了一周，终于算是看了个大概，知其然而不知其所以然...，与LR回归不一样的是这里的分类结果是多类，这里举一个网上很常见的例子，三种类型的点各100个，对其进行分类。
1、数据生成
①数据生成没什么难点，就是中间有点小细节需要处理
②scatter的参数都是不规则写法
③只有是一维的时候reshape和转置才相同，否则不一致（没有注意到，坑了很久）
def load_data():
    np.random.seed(0)    # 每次随机数固定
    N=100                #每类100个点
    D=2                  #每个点，(x,y)
    K=3                  #3类
    X=np.zeros((300,2))
    Y=np.zeros((300,1))
    for i in range(K):
        ix = range(N * i, N * (i + 1))   #0-100、100-200、200-300
        r = np.linspace(0.0, 1, N)  # radius
        t = np.linspace(i * 4, (i + 1) * 4, N) + np.random.randn(N) * 0.2
        X[ix] = np.c_[r * np.sin(t), r * np.cos(t)]
        Y[ix] = i
    #注意：这里的c=Y，Y要转换为不规则写法，s表示点的大小，cmap表示类型
    plt.scatter(X[:, 0], X[:, 1], c=Y.squeeze(), s=30, cmap=plt.cm.Spectral)
    #plt.scatter(X[:, 0], X[:, 1], c=Y.ravel(), s=30, cmap=plt.cm.Spectral)
    #plt.scatter(X[:, 0], X[:, 1], c=Y.reshape(300), s=30, cmap=plt.cm.Spectral)
    plt.show()
    #X.T和X.reshape(2,300)不一样
    #Y.T和Y.reshpe(1,300)是一样的
    #还一个是数据类型转换，Y在传播过程中要当下表使用,计算损失函数
    return D,K,X.T,Y.reshape(1,300).astype(np.int32)

2、没有隐藏层，直接就是正向传播
softmax的正向传播的输出也与LR回归不同，有多个输出，激活函数与其它的不同，损失函数也与其它不同。主要思想如下：
①计算出Z （这里的Z是个3*300的矩阵，每一列表示为每个点的输出）
②激活函数就是操作Z为e^Z，然后归一化（列归一），计算出每种分类的概率
③损失函数，根据标签Y，Y作下标，找到每一列对应位置的概率使其经可能最大，-log(P)就是使其经可能小（注意这里只取标签位置对应的概率，也就是取300个，计算-log，然后取平均值）
④为了防止过拟合，需要加上一个正则损失

Z=np.dot(W,X)+B
A=np.exp(Z)
A=A/np.sum(A,axis=0,keepdims=True)
L=-np.log(A[Y,range(M)])         #实际标签对应的值求log、
loss=np.sum(L)/M + 0.5*reg*np.sum(W*W)
3、反向传播
反向传播最难的就是求导，k表示标签对应的概率，链式求导，dL/dZk可以求出来，同理dL/dZi=Pi（k表示标签对应的下标，i表示非标签对应的下标）(L只和Pk有关，但是Pk和Zi有关，因为Σ中包含Zi)，综合可以求出dL/dZ


A[Y, range(M)] -=1   #A表示概率P，且此操作把对应标签的概率减1，其余概率不变
dw=np.dot(A,X.T)/M
dw+=reg*W
db=np.sum(A,axis=1,keepdims=True)/M
W -= rate * dw
B -= rate * db
4、源码（一个生成数据函数，一个传播函数，一个预测函数，一个效果展示画图函数）
import numpy as np
import matplotlib.pyplot as plt
#设置画图字体
from matplotlib.font_manager import FontProperties
def load_data():
    np.random.seed(0)    # 每次随机数固定
    N=100                #每类100个点
    D=2                  #每个点，(x,y)
    K=3                  #3类
    X=np.zeros((300,2))
    Y=np.zeros((300,1))
    for i in range(K):
        ix = range(N * i, N * (i + 1))   #0-100、100-200、200-300
        r = np.linspace(0.0, 1, N)  # radius
        t = np.linspace(i * 4, (i + 1) * 4, N) + np.random.randn(N) * 0.2
        X[ix] = np.c_[r * np.sin(t), r * np.cos(t)]
        Y[ix] = i
    print(X[:, 0].shape)
    #注意：这里的c=Y，Y要转换为不规则写法，s表示点的大小，cmap表示类型
    plt.scatter(X[:, 0], X[:, 1], c=Y.squeeze(), s=30, cmap=plt.cm.Spectral)
    #plt.scatter(X[:, 0], X[:, 1], c=Y.ravel(), s=30, cmap=plt.cm.Spectral)
    #plt.scatter(X[:, 0], X[:, 1], c=Y.reshape(300), s=30, cmap=plt.cm.Spectral)
    plt.show()
    #X.T和X.reshape(2,300)不一样
    #Y.T和Y.reshpe(1,300)是一样的
    #还一个是数据类型转换，Y在传播过程中要当下表使用,计算损失函数
    return D,K,X.T,Y.reshape(1,300).astype(np.int32)
def plot_answer(model,X,Y):
    x_min, x_max = X[0, :].min() - 1, X[0, :].max() + 1
    y_min, y_max = X[1, :].min() - 1, X[1, :].max() + 1
    h = 0.01
    # xx和yy是两个大小相等的矩阵
    # xx和yy是提供坐标（xx,yy）
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    # 列扩展，每个点是一行
    Z = model(np.c_[xx.ravel(),yy.ravel()].T)
    Z = Z.reshape(xx.shape)
    # 轮廓,等高线，预测值相等的点描绘成一个轮廓，把图像进行分割
    # contour   没有颜色填充，只是分割
    # contourf  会进行颜色填充
    plt.contourf(xx, yy, Z, cmap=plt.cm.Spectral)
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    plt.ylabel(u'y轴', fontproperties=font)
    plt.xlabel(u'x轴', fontproperties=font)
    plt.scatter(X[0, :], X[1, :], c=Y.squeeze(), s=30, cmap=plt.cm.Spectral)
    plt.show()
def propagation():
    M=X.shape[1]
    reg=0.001
    rate=1.0
    W=0.01*np.random.randn(K,D)
    B=np.zeros((K,1))
    for i in range(200):
        Z=np.dot(W,X)+B
        A=np.exp(Z)
        A=A/np.sum(A,axis=0,keepdims=True)
        L=-np.log(A[Y,range(M)])         #实际标签对应的值求log、
        loss=np.sum(L)/M + 0.5*reg*np.sum(W*W)
        if(i%10==0):
            print("iteration %d: loss %f" % (i, loss))
        A[Y, range(M)] -=1
        dw=np.dot(A,X.T)/M
        dw+=reg*W
        db=np.sum(A,axis=1,keepdims=True)/M
        W -= rate * dw
        B -= rate * db
    Z=np.dot(W,X)+B
    predict=np.argmax(Z,axis=0)
    print('training accuracy: %.0f%%' % (np.mean(predict == Y)*100))
    return W,B
def predict(X):
    Z=np.dot(W,X)+B
    ans=np.argmax(Z,axis=0)
    return ans
D,K,X,Y=load_data()
W,B=propagation()
plot_answer(lambda x: predict(x),X,Y)

总结：
①激活函数、损失函数、反向传播（激活函数就是先指数操作，再归一化求概率，损失函数就是-log，损失函数操作的是正确标签对应的概率，反向传播，链式求导，一步一步慢慢推）
②遇到的问题和学习到的知识（np.argmax、np.argsort，终于知道python的强大了，不需要二维排序，直接可以使用argmax求出最大概率的下标是多少、或者求出下标的排序结果。）（转置和reshape可不是同一个概念）（A[1][:]和A[1,:]的使用，遍历时支持后者）（标签Y的数据类型的转换，因为要当下标使用，A.dtype#查看数据类型，A为Numpy类型。A.astype(np.float64)#转换数据类型。A.astype(np.int32/np.int64)  #转换数据类型）
③准确率才49%，需要添加隐藏层才能进一步提高准确率

        差错控制技术主要是在编码端通过改变策略或添加冗余的方式，将错误限制在一定的范围之内防止误码扩散。使用该法当然是一定程度上保证了视频质量，却增加了带宽。网络流量很可能不够，速度会比较慢。相比之下，错误隐藏是一种在解码端进行处理的技术。错误隐藏就是解码器在不需要从编码器得到额外信息的前提下，利用当前帧已接收宏块或先前已接收视频序列的相关性对丢失宏块进行恢复，以主观可接受的近似原来质量的视频数据来隐藏错误受损的数据。错误隐藏技术要求算法尽可能简单，尤其是实时性比较高的场合。错误也分为帧间隐藏（intra
 frame，一帧的内部，spatial）和帧内隐藏（inter frame，空间temporal）。
        帧内错误比较容易处理，先来说一下帧内错误隐藏的算法实现：采用基于边缘检测错误隐藏或空间像素加权平均算法。


        先来说一下丢失块分成哪些类型：
        Lena：边缘信息比较丰富，局部纹理比较复杂；
        peppers：边缘信息丰富，纹理比较复杂；
        baboon：边缘信息较少，纹理复杂；
        forman：边缘信息丰富，纹理较少

     



        再来说一下边缘检测：
        Sobel边缘检测算子：
        http://baike.baidu.com/link?url=_tc38z77JqCqbe47WzeZxjsFvPnJJoeWGo2ZpHjkd78lfP1T54XDifUQ_NSeEuaplcyzajLrBq4vI-Zo8i6Jda
        http://www.cnblogs.com/lancidie/archive/2011/07/17/2108885.html
        http://homepages.inf.ed.ac.uk/rbf/HIPR2/featops.htm （这里包含了几乎边缘检测的算子，这个网站真是极好的）

        我们要知道怎么计算：
      
        我们可以把最中间的那个看做我们要求的点：

   该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像灰度值，其公式如下：
  
  
      Sobel算子根据像素点上下、左右邻点灰度加权差，在边缘处达到极值这一现象检测边缘。对噪声具有平滑作用，提供较为精确的边缘方向信息，边缘定位精度不够高。当对精度要求不是很高时，是一种较为常用的边缘检测方法。

   普利维特算子(Prewitt operate)：有噪声平滑作用，可是精度不高。
 

         
    罗伯茨交叉边缘检测（Roberts Cross operator）：精度高，可是计算量大，对噪声比较敏感。
   

            当然还有canny算子什么的，还想深入了解就看我给的那个地址了。


           实现基于边缘检测的单一方向的错误隐藏：
           （1）进行边缘检测，获取每个像素点的梯度及方向索引K（K由角度决定）。使用Sobel算子，对丢失块周围相邻像素进行检测，得到G和角度a，a取值为（-pi，pi）。对于角度小于0的要加上pi，再对pi/8进行模运算得到。
           （2）判断该点的边缘信息是否会通过丢失块，如果是则判断为相关像素，并将其梯度累加到该方向下的数组中去。要判断相关像素，通过角度θ及其坐标判断的。比如有两个点：p1（x1，y1），p2（x2，y2），方向都是θ，则它们在丢失块上边界的交点分别为X1 = 6，X2 = 17，计算公式为：
             X1 = xn + （8+1-yn）/ tan θ  像素检测向下延伸后能够通过下边界
             X2 = xn + （8+16-yn）/ tan θ 像素检测向上延伸后能够通过下边界

            Y1 = yn + （8+1-xn）/ tan θ  像素检测向右延伸后能够通过下边界
            Y2 = yn + （8+16-xn）/ tan θ 像素检测向左延伸后能够通过下边界
             这里的X1，X2，Y1，Y2分别代表丢失块在相邻像素中的边界上的某个像素点。X1的横坐标固定为8，X2的横坐标固定为25，Y1的纵坐标固定为8，Y2的纵坐标固定为25。当某个像素通过延伸后边界坐标在9~24时，可以判断该点为相关像素。就将该梯度累加到下标为SG[k]当中。如上图p2方向为K，SG[k] = SG[k] + G(xn,yn)
           （3）方向插值，根据最大梯度下的方向编号对丢失宏块进行像素插值。根据最大梯度对应一个K，用K作为丢失的插值方向Edgdir。
             插值要首先计算出选择相邻边界中那个像素点来进行插值，然后计算出插值点和这两个像素点之间的距离进行加权，最后完成插值算法。选择点的方法和像素判断方法一样。对丢失块中每个像素点进行K方向下延伸，判断其与X1，X2，Y1，Y2的交点，当其交点在边界范围之内时，则可以判断出用于插值的相邻边界像素点坐标。

             X1 = xn + （8+1-yn）/ tan （Edgdir *π/8）  像素检测向下延伸后能够通过下边界
             X2 = xn + （8+16-yn）/ tan（Edgdir *π/8） 像素检测向上延伸后能够通过下边界

            Y1 = yn + （8+1-xn）/ tan（Edgdir *π/8）像素检测向右延伸后能够通过下边界
            Y2 = yn + （8+16-xn）/ tan （Edgdir *π/8）像素检测向左延伸后能够通过下边界 

            (xn , yn)为丢失宏块中某个像素点的坐标位置；X1，X2，Y1，Y2为边界；Edgdir为方向编号。若得到的两个边界点在9~24范围内，取整之后就是可用于插值的相邻边界像素。
            比如我们得到两个点pp1，pp2，我们可利用公式：
            
            P，pp1，pp2为像素点的像素值；d1，d2是通过P，pp1，pp2三个像素点的坐标，通过两点距离公式计算得到的两点之间的距离。
        以上方法基于边缘检测错误隐藏算法只单独一个方向上插值的时候比较好！如果边缘信息太复杂的时候就会出现虚假边缘。这时候就采用空间算数加权算法可以得到比较好的效果。什么时候采用边缘检测的单方向什么时候采用空间像素加权平均算法。采用将最大梯度SG[k]比总共的梯度，我们设置一个阈值，比如0.2。若比值大于阈值，说明边缘信息比较简单，所以可以进行边缘检测的单方向线性插值。否则采用空间像素加权平均算法实现。根据这个比值，我们选择是采用空间像素加权平均还是边缘检测的单一方向插值算法。
        所以对于图像的边缘信息在一个方向下进行插值，因此对于边缘信息比较单一的丢失块隐藏效果较好。对于纹理复杂的丢失块，空间像素加权平均算法效果较好。




    H．264 
是ISO运动图像专家组与ITU-T的视频编码专家组联合推出的新一代视频压缩编码标准，其主要目的是提高压缩效率和增加网络友好性。由于使用了高效的编码技术，如帧内预测、1
 /4像素精度运动补偿、多参考帧等，H．264的
 编码效率得到极大提高，但也使 得编码后的码流对错误更加敏感，任何一个比特的错误不仅可能破坏当前帧，还可能造成错误扩散，导致后续视频帧质量的严重下降
。错误隐藏就是当图像发生错误时，在解码器端充分利用视频数据的相关性进行的图像修复措施
 ，以降低错误的影响，提高视频的主 观质量。H．264错误隐藏分为帧间错误隐藏和帧内错误隐藏，前者利用时间相关性
 ，通过先前已经正确接收到的图像帧来恢复当前丢失的图像帧信息。后者则是利用空间上的相关性
 ，通过一帧内已经接收到的宏块来恢复丢失宏块的信息。空间加权像素平均算法是H.264自带的帧内错误隐藏算法，该方法对平滑区域的隐藏效果较好，但是对于有边缘的信息块，隐藏效果较差。

参考论文：
基于H.264/AVC的自适应视频错误隐藏技术算法  包宋建等
基于边缘检测及方向加权的H.264帧内错误隐藏算法  李强等

【实验名称】基于龟壳的信息隐藏实验       

值得说明的是，网上有很多龟壳算法其实现过程是不够完整的，省去了一些例如边界处理的部分。而本报告中的代码对于龟壳算法实现的是完整的，可以说是一个比较完美的复现。

【实验目的】

1、加深对图像信息隐藏基本原理和方法的掌握

2、学习并实现基于龟壳的信息隐藏和提取算法

 

【实验原理】

1、龟壳矩阵为一种矩阵元素有着特殊分布规律的N*N（常见如256*256）矩阵，其元素值为0-7，按照一定的差值规律连续组成。因在矩阵中值为0-7的8个元素连线形状恰如龟壳而得名

2、基于龟壳魔术矩阵的方法修改载体图像的像素来隐藏秘密信息，每个像素平均可以隐藏1.5个比特。隐藏信息后的图像用人类的肉眼很难发现和原来的载体图像有任何差异

 

【实验内容】

代码：

function Turtle_Shell()
clc;clear;close all;
I=imread('coverlena.bmp');
figure;subplot(1,2,1);imshow(I);
title('原图');
Secret=23;                           %待隐藏数据的大小
data=round(rand(1,Secret));          %随机生成大小为Secret的待隐藏的2进制数据
fprintf('隐藏数据：');disp(data);
[length,width]=size(I);
data1=zeros(1,fix(Secret/3)+1);
point=size(data1,2);
for i=Secret:-3:1
    if i-3<0
        num=1;
        for i=i:-1:1
        data1(1,point)=data1(1,point)+num*data(1,i);
        num=2*num;
        end
        break;
    end  
    data1(1,point)= 4*data(1,i-2)+2*data(1,i-1)+data(1,i);
    point=point-1;
end
fprintf('转换成8进制：');disp(data1);
%生成龟壳矩阵
tortoise=Tortoise();

%进行信息隐藏
point=1;
for i=1:length
    if point>size(data1,2)
            break;
    end
    for j=1:2:width
        if point>size(data1,2)
            break;
        end
        %像素值0到255龟壳矩阵1到256，需将两者对齐,故此处tortoise输出时，行列+1
        %输出原始信息
        fprintf('待隐藏信息：%d  原像素值：%3d %3d  对应龟壳矩阵值：%d  ',...
            data1(point),I(i,j),I(i,j+1),tortoise(double(I(i,j))+1,double(I(i,j+1))+1));
        %进行信息隐藏，即像素值修改
        [I(i,j),I(i,j+1)]=find(I(i,j),I(i,j+1),data1(point),tortoise);
        %输出信息隐藏后信息
        fprintf('信息隐藏后像素值：%3d %3d  对应龟壳矩阵值：%d\n',...
            I(i,j),I(i,j+1),tortoise(double(I(i,j))+1,double(I(i,j+1))+1));
        
        point=point+1;  
    end
end
fprintf('信息隐藏成功!\r\n');
subplot(1,2,2);imshow(I);
title('隐藏后图像');

%隐藏信息读取
point=1;
for i=1:length
    if point>size(data1,2)
            break;
    end
    for j=1:2:width
        if point>size(data1,2)
            break;
        end
        %在龟壳矩阵中找出信息隐藏点，注意到像素值0到255，龟壳矩阵1到256，需将两者对齐
        Data1(point)=tortoise(double(I(i,j))+1,double(I(i,j+1))+1);
        point=point+1;  
    end
end

Data1=double(Data1);
fprintf('提取数据：');disp(Data1);

point=size(Data1,2);
for i=Secret:-3:1
    if i-3<0
        for i=i:-1:1
            Data(1,i)=mod(Data1(1,point),2);
            Data1(1,point)=(Data1(1,point)-mod(Data1(1,point),2))/2;
        end
        break;
    end  
    Data(1,i)=mod(Data1(1,point),2);
    Data1(1,point)=(Data1(1,point)-mod(Data1(1,point),2))/2;
    Data(1,i-1)=mod(Data1(1,point),2);
    Data1(1,point)=(Data1(1,point)-mod(Data1(1,point),2))/2;
    Data(1,i-2)=mod(Data1(1,point),2);
    point=point-1;
end
fprintf('转换成2进制：');disp(Data);
end

%生成龟壳矩阵
function tortoise=Tortoise()     
tortoise(1,1)=6;
for i=2:256
    tortoise(1,i)=mod(tortoise(1,i-1)+1,8);
end
for i=2:256
    if mod(i,2)==0
        tortoise(i,:)=mod(tortoise(i-1,:)-2,8);
    else
        tortoise(i,:)=mod(tortoise(i-1,:)-3,8);
    end
end
end

%进行信息隐藏,带入当前两个像素值，返回信息隐藏点的行列
function [X,Y]=find(X,Y,number,tortoise)   
%像素值0到255，龟壳矩阵1到256，需将两者对齐
X=double(X);Y=double(Y);
X=X+1;Y=Y+1;

%对边界进行特殊处理
if (X==1&&Y==1)||(X==256&&Y==1)||(X==1&&Y==256)||(X==256&&Y==256)||...
        (X==1&&mod(Y,2)==1)||(X==256&&mod(Y,2)==1)||...
        (Y==1&&(mod(X,4)==0||mod(X,4)==1))||(Y==256&&(mod(X,4)==2||mod(X,4)==3))
    fprintf('该点在边界上    ');
    if X==1&&Y==1             %判断其是否是龟壳矩阵的左上顶点
        [X1,Y1,X2,Y2]=matrix(X,Y,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif X==256&&Y==1       %判断其是否是龟壳矩阵的左下顶点
        [X1,Y1,X2,Y2]=matrix(X-2,Y,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif X==1&&Y==256       %判断其是否是龟壳矩阵的右上顶点
        [X1,Y1,X2,Y2]=matrix(X,Y-2,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif X==256&&Y==256     %判断其是否是龟壳矩阵的右下顶点
        [X1,Y1,X2,Y2]=matrix(X-2,Y-2,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif X==1&&mod(Y,2)==1  %判断其是否在龟壳矩阵上边界、是否需要边界处理
        [X1,Y1,X2,Y2]=matrix(X,Y-1,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif X==256&&mod(Y,2)==1%判断其是否在龟壳矩阵下边界、是否需要边界处理
        [X1,Y1,X2,Y2]=matrix(X-2,Y-1,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif Y==1&&(mod(X,4)==0||mod(X,4)==1)%判断其是否在龟壳矩阵左边界、是否需要边界处理
        [X1,Y1,X2,Y2]=matrix(X-1,Y,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    elseif Y==256&&(mod(X,4)==2||mod(X,4)==3)%判断其是否在龟壳矩阵右边界、是否需要边界处理
        [X1,Y1,X2,Y2]=matrix(X-1,Y-2,number,tortoise);
        temp1=X1;temp2=Y1;
        if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
            temp1=X2;temp2=Y2;
        end
        X=temp1;Y=temp2;
    end
    X=X-1;Y=Y-1;  %将像素值恢复到0到255范围内
    return;
end

%判断其是否是龟壳内部两点中的上内部点
if (tortoise(X,Y)==1 && tortoise(X+1,Y)==6)||(tortoise(X,Y)==3 && tortoise(X+1,Y)==0)||...
        (tortoise(X,Y)==5 && tortoise(X+1,Y)==2)||(tortoise(X,Y)==7 && tortoise(X+1,Y)==4)
    fprintf('该点在龟壳内    ');
    [X,Y]=traversal(X-1,Y,number,tortoise);
    X=X-1;Y=Y-1;  %将像素值恢复到0到255范围内
    return;
end

%判断其是否是龟壳内部两点中的下内部点
if (tortoise(X,Y)==6 && tortoise(X-1,Y)==1)||(tortoise(X,Y)==0 && tortoise(X-1,Y)==3)||...
        (tortoise(X,Y)==2 && tortoise(X-1,Y)==5)||(tortoise(X,Y)==4 && tortoise(X-1,Y)==7)
    fprintf('该点在龟壳内    ');
    [X,Y]=traversal(X-2,Y,number,tortoise);
    X=X-1;Y=Y-1;  %将像素值恢复到0到255范围内
    return;
end

%剩下的点全在龟壳的边上，且只有上下顶点两类
if mod(tortoise(X,Y),2)==1          %判断其是否在龟壳边的上顶点
    fprintf('该点在龟壳边    ');
    [X1,Y1]=traversal(X,Y,number,tortoise);
    [X2,Y2]=traversal(X-2,Y-1,number,tortoise);
    [X3,Y3]=traversal(X-2,Y+1,number,tortoise);
    temp1=X1;temp2=Y1;
    if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
        temp1=X2;temp2=Y2;
    end
    if (X3-X)^2+(Y3-Y)^2<(temp1-X)^2+(temp2-Y)^2
        temp1=X3;temp2=Y3;
    end
    X=temp1;Y=temp2;
    X=X-1;Y=Y-1;  %将像素值恢复到0到255范围内
    return;
end

if mod(tortoise(X,Y),2)==0          %判断其是否在龟壳边的下顶点
    fprintf('该点在龟壳边    ');
    [X1,Y1]=traversal(X-3,Y,number,tortoise);
    [X2,Y2]=traversal(X-1,Y-1,number,tortoise);
    [X3,Y3]=traversal(X-1,Y+1,number,tortoise);
    temp1=X1;temp2=Y1;
    if (X2-X)^2+(Y2-Y)^2<(temp1-X)^2+(temp2-Y)^2
        temp1=X2;temp2=Y2;
    end
    if (X3-X)^2+(Y3-Y)^2<(temp1-X)^2+(temp2-Y)^2
        temp1=X3;temp2=Y3;
    end
    X=temp1;Y=temp2;
    X=X-1;Y=Y-1;  %将像素值恢复到0到255范围内
    return;
end
end

%从上顶点开始遍历一个龟壳，并返回需要寻找点的坐标
function [X,Y]=traversal(X,Y,number,tortoise)
if X<1||X>256||X+1<1||X+1>256||X+2<1||X+2>256||X+3<1||X+3>256||...
        Y<1||Y>256||Y-1<1||Y-1>256||Y+1<1||Y+1>256
    X=inf;Y=inf;
    elseif tortoise(X,Y)==number
    elseif tortoise(X+1,Y)==number
        X=X+1;
    elseif tortoise(X+1,Y-1)==number
        X=X+1;
        Y=Y-1;
    elseif tortoise(X+1,Y+1)==number
        X=X+1;
        Y=Y+1;
    elseif tortoise(X+2,Y)==number
        X=X+2;
    elseif tortoise(X+2,Y-1)==number
        X=X+2;
        Y=Y-1;
    elseif tortoise(X+2,Y+1)==number
        X=X+2;
        Y=Y+1;
    elseif tortoise(X+3,Y)==number
        X=X+3;
end
end

%分别从左上顶点和右下顶点开始，遍历一个3*3矩阵，并返回需要寻找点的坐标
%分别从左上顶点和右下顶点开始遍历的作用是找出距离最小的点
function [X1,Y1,X2,Y2]=matrix(X1,Y1,number,tortoise)
flag=0;
X2=X1;Y2=Y1;
for i=0:1:2
    for j=0:1:2
        if tortoise(X1+i,Y1+j)==number
            X1=X1+i;
            Y1=Y1+j;
            flag=1;
            break;
        end
    end
    if flag==1
        break;
    end
end
flag=0;
for i=2:-1:0
    for j=2:-1:0
        if tortoise(X2+i,Y2+j)==number
            X2=X2+i;
            Y2=Y2+j;
            flag=1;
            break;
        end
    end
    if flag==1
        break;
    end
end
end


实验结果：





 

其中，生成的龟壳矩阵（部分）如下：



分析实验代码和结果，可以看到：

整个过程是先随机生成10个2进制数，并将其转换成8进制数，再进行信息隐藏，最后进行隐藏数据的提取，观察是否隐藏成功。

在该示例结果中，既有点在龟壳内部（第4个待隐藏值）也有点在龟壳边上（第1、2、3个待隐藏值），既有点对应的像素值不是待隐藏值（第1、3、4个待隐藏值）也有点刚好是待隐藏值（第2个待隐藏值）。实验结果显示，待隐藏信息全都顺利隐藏并且成功读取。

而且根据龟壳矩阵进行手动分析，可以验证整个实验过程是正确的。

 

但是由于图像“coverlena.bmp”的特殊性，实验中并不会遇到点在边界上的情况。所以，我手动构造了一个特殊的8×8图像，以检查其边界的特判处理是否正确。

这个图像包括：龟壳矩阵四个顶点（图像第1行）、龟壳矩阵上下左右四条边上随机各取一点（图像第2行）。图像具体像素值如下：



 

实验结果：

 





注意到其中第8、9个待隐藏信息显示在龟壳边上而不在边界上，这是因为虽然取的这八个点都是在龟壳矩阵的边上，却并不意味着需要进行“特判处理”即不在“边界”上，龟壳矩阵的四条边上也有一部分点是在龟壳上的。

实验结果显示，待隐藏信息全都顺利隐藏并且成功读取。根据龟壳矩阵进行手动分析，可以验证整个实验过程是正确的。

至此，对龟壳内、龟壳边和边界上的点都进行了测试，实验结果和手动验证都显示算法是正确的。

 

【小结或讨论】

基于龟壳的信息隐藏算法从本质上来说和EMD算法是一样的，都是利用魔术矩阵进行信息隐藏。但EMD的有一个很大的缺点，就是它隐藏的数据需要是５进制的，这无疑会使信息隐藏的效率变得低下。而基于龟壳的信息隐藏算法很好地解决了这个问题，它采用的是８进制的魔术矩阵，２进制同８进制的转换远比同５进制的转换有效率的多，而且使得每个信息隐藏点的容量进一步扩大到3bits。

基于龟壳的信息隐藏算法的基本原理：

龟壳矩阵是这个算法的基础，其中每个元素的值都在0-7之间，并且按照一定的差值规律排列，每一个龟壳都无冗余的包含了0-7的8个数。在进行信息隐藏时，先以图像中的两个像素值作为横、纵坐标来找到龟壳矩阵中对应的元素值，再在对应的龟壳中找到待隐藏信息值相同的点，以该点的横纵坐标替换掉数字图像中的原像素值。这样使得隐藏容量增加的同时，图像失真却保持很小。

当然效率的大幅提升必然会导致算法复杂度的提升。基于龟壳的信息隐藏算法和上一次实验的LSB和EMD算法相比难度大大增加，导致整个算法完整实现后，代码近300行。由于算法的复杂，我想在这里总结一下整个代码的构写思路和最核心的在龟壳中寻找待隐藏值的思路。

代码构写思路：

整个代码包含以下5个函数

1、Turtle_Shell()为主函数

2、Tortoise()作用是生成龟壳矩阵

3、find(X,Y,number,tortoise)作用是找到对应龟壳并找到待隐藏点，进行像素值替换。其中X为行数，Y为列数，number为待隐藏值，tortoise是已生成的龟壳矩阵。之所以将龟壳矩阵带入，是为了避免每次进行信息隐藏时都重新生成龟壳矩阵，从而降低代码运行效率。

4、traversal(X,Y,number,tortoise)作用是从上顶点开始遍历一个龟壳。其中X为该龟壳上顶点行数，Y为该龟壳上顶点列数。比较重要的一点是，若遍历龟壳不完整则返回无穷大，这是为了靠近边界处有时只需要遍历1个或2个龟壳而设立的。其他同上。

5、matrix(X1,Y1,number,tortoise)作用是分别从矩阵的左上顶点和右下顶点开始，遍历一个3×3矩阵，并返回需要寻找点的坐标。分别从两个起点的作用是两者遍历的路径刚好互补，从而若与待隐藏值相同的点有两个，可以找出距离最小的点。其他同上。

主函数Turtle_Shell控制整个流程，先调用Tortoise生成龟壳矩阵，再通过find函数进行信息隐藏。在find函数中，分类判断了待隐藏值的位置后，若在龟壳边/内则调用traversal遍历周围对应龟壳，若在边界则调用matrix遍历3×3矩阵，最终找到信息隐藏点并返回行列值。

在龟壳中寻找待隐藏值的思路：

1、判断点在边界的特殊情况

（1）、四个顶点

（2）、龟壳矩阵第1行且列数为奇数（上边界特判）

（2）、龟壳矩阵第256行且列数为奇数（下边界特判）

（2）、龟壳矩阵第1列且行数为模4余0或1（左边界特判）

（2）、龟壳矩阵第256列且行数为模4余2或3（右边界特判）

2、判断点在龟壳内的情况

根据观察，整个魔术矩阵有4种龟壳，其内部两个值各不相同，分别为：<1,6>、<3,0>、<5,2>、<7,4>。若符合其中一种排列情况，说明其在龟壳内部。

3、判断点在龟壳边的情况

排除了以上两种情况后，该点则必在龟壳边上。根据分析，在龟壳边上的点实际上只有两类：龟壳上顶点和龟壳下顶点。

（1）、该点值为奇数则为龟壳上顶点

（2）、该点值为偶数则为龟壳下顶点

这次实验比较复杂，整个算法也是在不断地试错和分析中逐步完善的。实验过程中遇到过很多bug，比如我在将像素值和龟壳矩阵对齐时，遇见255+1始终等于255的奇怪现象。最后分析才发现，原来像素值默认为uint8格式，其表示的最大值就是255，最后将变量转换成double类型就解决了问题。