Somogyi Stucki dithering的matalab代码

以下是Somogyi Stucki抖动算法的MATLAB代码：

function output_img = dither_somogyi_stucki(input_img)
% Somogyi Stucki dithering algorithm
% input_img: grayscale input image (double, range 0-1)
% output_img: dithered image (binary)
    
    % Define the dithering matrix
    dither_mat = [0 0 0 8 4; 2 4 8 4 2; 1 2 4 2 1]/42;
    
    % Define the dimensions of the input image
    [height, width] = size(input_img);
    
    % Initialize the output image
    output_img = zeros(height, width);
    
    % Iterate over each pixel
    for y = 1:height
        for x = 1:width
            % Calculate the current pixel's value
            old_pixel = input_img(y,x);
            
            % Calculate the new pixel's value
            new_pixel = old_pixel + dither_mat(mod(y-1,3)+1,mod(x-1,5)+1);
            
            % Threshold the new pixel's value
            if new_pixel >= 0.5
                output_img(y,x) = 1;
            else
                output_img(y,x) = 0;
            end
            
            % Calculate the quantization error
            quant_error = old_pixel - output_img(y,x);
            
            % Distribute the quantization error to neighboring pixels
            if x < width
                input_img(y,x+1) = input_img(y,x+1) + quant_error * 8/42;
            end
            if x < width-1
                input_img(y,x+2) = input_img(y,x+2) + quant_error * 4/42;
            end
            if x > 1 && y < height
                input_img(y+1,x-1) = input_img(y+1,x-1) + quant_error * 2/42;
            end
            if y < height
                input_img(y+1,x) = input_img(y+1,x) + quant_error * 4/42;
            end
            if y < height-1
                input_img(y+2,x) = input_img(y+2,x) + quant_error * 2/42;
            end
            if x < width && y < height
                input_img(y+1,x+1) = input_img(y+1,x+1) + quant_error * 8/42;
            end
            if x < width && y < height-1
                input_img(y+2,x+1) = input_img(y+2,x+1) + quant_error * 4/42;
            end
            if x < width-1 && y < height-1
                input_img(y+2,x+2) = input_img(y+2,x+2) + quant_error * 2/42;
            end
        end
    end
    
    % Convert the output image to binary
    output_img = logical(output_img);
end

请注意，该代码假定输入图像是灰度图像，并且其值范围为0到1。如果你的输入图像不是灰度图像或其值范围不是0到1，你需要先进行适当的处理。