msrcr算法matlab代码

下面是msrcr算法的MATLAB代码：

function [I_MS_out, I_PAN_out] = msrcr(I_MS, I_PAN, ratio)
% Multi-Scale Retinex with Color Restoration (MSRCR)
% Inputs:
%   I_MS:   multi-spectral image, MxNxL
%   I_PAN:  panchromatic image, MxN
%   ratio:  ratio of the weight between the MSRCR output and the input MS image
% Outputs:
%   I_MS_out:   enhanced MS image, MxNxL
%   I_PAN_out:  downsampled and enhanced PAN image, MxN
 
% The code is modified based on the following paper:
%   Rahman Z, Jobson DJ, Woodell GA. Multi-scale retinex for color image 
%   enhancement. In Proceedings of the International Conference on Image 
%   Processing 1996 Oct 27 (Vol. 3, pp. 1003-1006). IEEE.
%--------------------------------------------------------------------------
% Initial settings
if nargin == 2
    ratio = 0.5;
end

[M,N,L] = size(I_MS);
  
% Apply MSRCR to each spectral band
I_MS_out = zeros(size(I_MS));
for i=1:L
    I_MS_out(:,:,i) = MSRCR(I_MS(:,:,i));
end

% Apply MSRCR to the panchromatic band
I_PAN_out = MSRCR(I_PAN);

% Color restoration
I_MS_out = (1-ratio)*I_MS + ratio*I_MS_out;
I_PAN_out = (1-ratio)*I_PAN + ratio*I_PAN_out;
I_PAN_out = imresize(I_PAN_out,[M,N]);

%--------------------------------------------------------------------------
% MSRCR function
function I_out = MSRCR(I)

% Parameters
sigma_list = [15,80,250];
G = fspecial('gaussian',[15 15],7);

% Multi-scale processing
I_log = log(double(I)+1);
I_out = zeros(size(I));
for i=1:length(sigma_list)
    % Gaussian filtering
    I_smooth = imfilter(I_log,G*sigma_list(i),'same','replicate');
    % Subtract the smoothed version from the original image
    I_detail = I_log - I_smooth;
    % Rescale the detail so that it has a mean of 0.5 and a standard deviation of 0.25
    I_detail = (I_detail - mean2(I_detail)) / std2(I_detail) * 0.25 + 0.5;
    % Add the rescaled detail back to the smoothed image
    I_out = I_out + I_detail;
end
I_out = exp(I_out - log(length(sigma_list)+eps));

% Histogram equalization
I_out = histeq(I_out,1024);

%--------------------------------------------------------------------------

使用方法：

% 读取MS和PAN图像
I_MS = imread('ms_image.tif');
I_PAN = imread('pan_image.tif');

% 调用MSRCR函数进行图像增强
[I_MS_out, I_PAN_out] = msrcr(I_MS, I_PAN, 0.5);

% 显示增强后的MS和PAN图像
figure;imshow(I_MS_out);
figure;imshow(I_PAN_out);