有界广义高斯混合模型的图像去噪代码 matlab

以下是基于有界广义高斯混合模型的图像去噪的 MATLAB 代码：

function [denoised_img] = bggmm_denoising(noisy_img, sigma, k, alpha)

% 读入图像并归一化
img = im2double(noisy_img);

% 初始化参数
[height, width] = size(img);
patch_size = 8;
step_size = 4;
num_patches = floor((height - patch_size) / step_size + 1) * floor((width - patch_size) / step_size + 1);
patches = zeros(patch_size * patch_size, num_patches);

% 提取图像的重叠补丁
patch_index = 1;
for i = 1 : step_size : height - patch_size + 1
    for j = 1 : step_size : width - patch_size + 1
        patch = img(i : i + patch_size - 1, j : j + patch_size - 1);
        patches(:, patch_index) = patch(:);
        patch_index = patch_index + 1;
    end
end

% 用BG-GMM模型去噪
denoised_patches = zeros(patch_size * patch_size, num_patches);
for i = 1 : num_patches
    patch = patches(:, i);

    % 初始化BG-GMM模型
    bggmm = bggmm_init(patch, k, alpha);

    % 估计BG-GMM模型参数
    [bggmm, ~] = bggmm_estimate(patch, bggmm, sigma, 20);

    % 计算后验概率
    bg_prob = bggmm.w' * bggmm.p;
    fg_prob = 1 - bg_prob;
    bg_prob = reshape(bg_prob, [patch_size, patch_size]);
    fg_prob = reshape(fg_prob, [patch_size, patch_size]);

    % 用后验概率加权的均值来计算去噪后的补丁
    denoised_patch = sum(bsxfun(@times, bggmm.mu, bg_prob(:)), 2);
    denoised_patches(:, i) = denoised_patch;
end

% 重构图像
denoised_img = zeros(height, width);
count = zeros(height, width);
patch_index = 1;
for i = 1 : step_size : height - patch_size + 1
    for j = 1 : step_size : width - patch_size + 1
        denoised_patch = reshape(denoised_patches(:, patch_index), [patch_size, patch_size]);
        denoised_img(i : i + patch_size - 1, j : j + patch_size - 1) = denoised_img(i : i + patch_size - 1, j : j + patch_size - 1) + denoised_patch;
        count(i : i + patch_size - 1, j : j + patch_size - 1) = count(i : i + patch_size - 1, j : j + patch_size - 1) + 1;
        patch_index = patch_index + 1;
    end
end
denoised_img = denoised_img ./ count;
denoised_img = im2uint8(denoised_img);

end

function [bggmm] = bggmm_init(patch, k, alpha)

% 初始化BG-GMM模型
bggmm.k = k;
bggmm.alpha = alpha;
bggmm.mu = zeros(size(patch, 1), k);
bggmm.sigma = repmat(eye(size(patch, 1)) * sqrt(0.1), [1, 1, k]);
bggmm.w = ones(1, k) / k;
bggmm.p = ones(size(patch, 1), k) / k;

% 随机初始化均值
for i = 1 : k
    bggmm.mu(:, i) = patch(:, randi(size(patch, 2)));
end

end

function [bggmm, llh] = bggmm_estimate(patch, bggmm, sigma, max_iter)

% 估计BG-GMM模型参数
for iter = 1 : max_iter
    % 计算后验概率
    bg_prob = bggmm.w' * bggmm.p;
    fg_prob = 1 - bg_prob;
    bg_prob(bg_prob == 0) = 1e-10;
    fg_prob(fg_prob == 0) = 1e-10;
    bg_prob = reshape(bg_prob, [size(patch, 1), size(patch, 2)]);
    fg_prob = reshape(fg_prob, [size(patch, 1), size(patch, 2)]);

    % 更新均值和协方差矩阵
    for i = 1 : bggmm.k
        % 计算加权样本点数
        N = sum(bg_prob(:));
        if N < eps
            break;
        end

        % 更新均值
        bggmm.mu(:, i) = sum(bsxfun(@times, patch, bg_prob(:, i)), 2) ./ N;

        % 更新协方差矩阵
        diff = bsxfun(@minus, patch, bggmm.mu(:, i));
        bggmm.sigma(:, :, i) = sum(bsxfun(@times, diff, permute(diff, [1, 3, 2])) .* bg_prob(:, :, ones(size(patch, 1), 1)), 3) ./ N;

        % 防止协方差矩阵奇异
        bggmm.sigma(:, :, i) = bggmm.sigma(:, :, i) + eye(size(patch, 1)) * sigma^2;
    end

    % 更新后验概率和混合系数
    for i = 1 : bggmm.k
        % 计算样本点权重
        bggmm.p(:, i) = bg_prob(:, i) ./ (bg_prob(:, i) + fg_prob(:, i) * bggmm.alpha);

        % 更新混合系数
        bggmm.w(i) = sum(bg_prob(:, i)) / numel(patch);
    end

    % 计算对数似然
    llh(iter) = sum(log(sum(bsxfun(@times, bggmm.w, bg_prob), 2)));

    % 如果对数似然收敛，则停止迭代
    if iter > 1 && llh(iter) - llh(iter - 1) < eps
        break;
    end
end

end

其中，`noisy_img` 是需要去噪的图像，`sigma` 是噪声标准差，`k` 是BG-GMM模型中高斯分量的数量，`alpha` 是BG-GMM模型中前景和背景差异的参数。函数返回去噪后的图像 `denoised_img`。

使用方法：

noisy_img = imread('noisy_img.png'); % 读入图像
sigma = 20; % 噪声标准差
k = 3; % BG-GMM高斯分量数量
alpha = 2.5; % 前景和背景差异参数
denoised_img = bggmm_denoising(noisy_img, sigma, k, alpha); % 图像去噪
imshow(denoised_img); % 显示去噪后的图像