原始载体图像为灰度图像，水印为黑白二值小图像，写出嵌入水印和提取水印步骤， Matlab编程实现算法。同屏显示原始图像、含水印图像和嵌入的水印，计算含水印图像的 PSNR； 同屏显示原始水印和提取的水印图像，计算提取水印的 NC 和 BER。

嵌入水印步骤：

1. 读取原始载体图像和水印图像，并将水印图像进行二值化处理。
2. 将原始载体图像分成若干个非重叠的小块，对每个小块进行离散余弦变换（DCT）。
3. 对每个小块的DCT系数进行量化，并根据水印图像的像素值将某些DCT系数的量化值加上或减去一个固定的值。
4. 将加密后的DCT系数通过逆DCT变换得到含水印图像。

提取水印步骤：

1. 读取含水印图像，并将其分成若干个非重叠的小块。
2. 对每个小块进行DCT变换，并提取出被加密的DCT系数。
3. 根据加密规则反推出加密前的DCT系数，并计算与加密前的量化值的差值。
4. 根据差值的符号判断水印像素值的0或1，并将其组合得到水印图像。

Matlab编程实现：

嵌入水印：

% 读取原始载体图像和水印图像，并将水印图像进行二值化处理
img = imread('original_img.bmp');
watermark = imread('watermark.bmp');
watermark = im2bw(watermark);

% 定义DCT变换的块大小和量化矩阵
block_size = 8;
quant_matrix = [16 11 10 16 24 40 51 61;
               12 12 14 19 26 58 60 55;
               14 13 16 24 40 57 69 56;
               14 17 22 29 51 87 80 62;
               18 22 37 56 68 109 103 77;
               24 35 55 64 81 104 113 92;
               49 64 78 87 103 121 120 101;
               72 92 95 98 112 100 103 99];

% 将原始载体图像分成若干个非重叠的小块，对每个小块进行离散余弦变换（DCT）
[height, width] = size(img);
watermarked_img = zeros(height, width);
for i = 1:block_size:height
    for j = 1:block_size:width
        % 取出一个小块
        block = img(i:i+block_size-1, j:j+block_size-1);
        % 进行DCT变换
        dct_block = dct2(block);
        % 对DCT系数进行量化
        quantized_block = round(dct_block./quant_matrix);
        % 根据水印图像的像素值将某些DCT系数的量化值加上或减去一个固定的值
        for k = 1:block_size
            for l = 1:block_size
                if watermark(i+k-1, j+l-1) == 0
                    if mod(k+l, 2) == 0
                        quantized_block(k, l) = quantized_block(k, l) - 1;
                    else
                        quantized_block(k, l) = quantized_block(k, l) + 1;
                    end
                else
                    if mod(k+l, 2) == 0
                        quantized_block(k, l) = quantized_block(k, l) + 1;
                    else
                        quantized_block(k, l) = quantized_block(k, l) - 1;
                    end
                end
            end
        end
        % 将加密后的DCT系数通过逆DCT变换得到含水印图像
        watermarked_block = uint8(idct2(quantized_block.*quant_matrix));
        watermarked_img(i:i+block_size-1, j:j+block_size-1) = watermarked_block;
    end
end

% 同屏显示原始图像、含水印图像和嵌入的水印，计算含水印图像的PSNR
original_psnr = psnr(img, watermarked_img);
figure;
subplot(1,3,1);
imshow(img);
title('Original Image');
subplot(1,3,2);
imshow(watermarked_img);
title('Watermarked Image');
subplot(1,3,3);
imshow(watermark);
title('Embedded Watermark');
watermarked_psnr = psnr(img, watermarked_img);
disp(['Original PSNR: ', num2str(original_psnr), ', Watermarked PSNR: ', num2str(watermarked_psnr)]);

提取水印：

% 读取含水印图像
watermarked_img = imread('watermarked_img.bmp');

% 将含水印图像分成若干个非重叠的小块
[height, width] = size(watermarked_img);
watermark_extracted = zeros(height, width);
for i = 1:block_size:height
    for j = 1:block_size:width
        % 取出一个小块
        block = watermarked_img(i:i+block_size-1, j:j+block_size-1);
        % 进行DCT变换
        dct_block = dct2(block);
        % 提取出被加密的DCT系数
        encrypted_dct = zeros(block_size);
        for k = 1:block_size
            for l = 1:block_size
                if mod(k+l, 2) == 0
                    if dct_block(k, l) > dct_block(k, l+1)
                        encrypted_dct(k, l) = 0;
                    else
                        encrypted_dct(k, l) = 1;
                    end
                else
                    if dct_block(k, l) > dct_block(k+1, l)
                        encrypted_dct(k, l) = 0;
                    else
                        encrypted_dct(k, l) = 1;
                    end
                end
            end
        end
        % 根据加密规则反推出加密前的DCT系数，并计算与加密前的量化值的差值
        extracted_block = quant_matrix.*encrypted_dct;
        % 进行逆DCT变换
        watermark_extracted(i:i+block_size-1, j:j+block_size-1) = uint8(idct2(extracted_block));
    end
end

% 同屏显示原始水印和提取的水印图像，计算提取水印的NC和BER
original_watermark = imread('watermark.bmp');
figure;
subplot(1,2,1);
imshow(original_watermark);
title('Original Watermark');
subplot(1,2,2);
imshow(watermark_extracted);
title('Extracted Watermark');
nc = corr2(original_watermark, watermark_extracted);
ber = sum(xor(original_watermark(:), watermark_extracted(:)))/numel(original_watermark);
disp(['NC: ', num2str(nc), ', BER: ', num2str(ber)]);