![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021695/image.png)
Problem in getting proper output image format
2 views (last 30 days)
Show older comments
Here in this, I used the script to apply Multi-Scale Retinex(MSR) to L channel of the image which is obtained from RGB to CIELAB converted. After applying MSR, I could not view the image as expected. Expecting suggestions in this.
% Multi-Scale Retinex on L Component of CIE-LAB Space
close all;
clc;
%I=imread("D:\Search Works\Underwater Dataset\raw-890 images\23_img_.png");
I=imread("23_img_.png");
Ilab = rgb2lab(I);
Ir = Ilab(:,:,1);
Ig = Ilab(:,:,2);
Ib = Ilab(:,:,3);
%%
%%%%%%%%%%Set the required parameters%%%%%%
G = 192;
b = -30;
alpha = 125;
beta = 46;
Ir_double=double(Ir);
%%
%%%%%%%%%%Set the Gaussian parameter%%%%%%
sigma_1=15; %Three Gaussian Kernels
sigma_2=80;
sigma_3=250;
[x, y]=meshgrid((-(size(Ir,2)-1)/2):(size(Ir,2)/2),(-(size(Ir,1)-1)/2):(size(Ir,1)/2));
gauss_1=exp(-(x.^2+y.^2)/(2*sigma_1*sigma_1)); %Calculate the Gaussian function
Gauss_1=gauss_1/sum(gauss_1(:)); %Normalization
gauss_2=exp(-(x.^2+y.^2)/(2*sigma_2*sigma_2));
Gauss_2=gauss_2/sum(gauss_2(:));
gauss_3=exp(-(x.^2+y.^2)/(2*sigma_3*sigma_3));
Gauss_3=gauss_3/sum(gauss_3(:));
%%
%%%%%%%%%%Operates on R component%%%%%%%
% MSR Section
Ir_log=log(Ir_double+1); %Converts an image to a logarithm field
f_Ir=fft2(Ir_double); %The image is Fourier transformed and converted to the frequency domain
%sigma = 15 processing results
fgauss=fft2(Gauss_1,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss); %Move the center of the frequency domain to zero
Rr=ifft2(fgauss.*f_Ir); %After convoluting, transform back into the airspace
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr1=Ir_log-Rr_log;
%sigma=80
fgauss=fft2(Gauss_2,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr2=Ir_log-Rr_log;
%sigma=250
fgauss=fft2(Gauss_3,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr3=Ir_log-Rr_log;
Rr=0.33*Rr1+0.34*Rr2+0.33*Rr3; %Weighted summation
MSR1 = Rr;
SSR1 = Rr2;
%SSR
min1 = min(min(SSR1));
max1 = max(max(SSR1));
SSR1 = uint8(255*(SSR1-min1)/(max1-min1));
%MSR
min1 = min(min(MSR1));
max1 = max(max(MSR1));
MSR1 = uint8(255*(MSR1-min1)/(max1-min1));
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ssr = cat(3,SSR1,Ig,Ib);
ssr = double(ssr);
ssr = lab2rgb(ssr);
msr = cat(3,MSR1,Ig,Ib);
msr = double(msr);
msr = lab2rgb(msr);
subplot(3,2,1);imshow(I);title('Original') %Show the original image
subplot(3,2,2);imshow(ssr);title('SSR')
subplot(3,2,3);imshow(msr);title('MSR')
0 Comments
Accepted Answer
Image Analyst
on 4 Jun 2022
Try this:
% Initialization steps.
clc;
clearvars;
close all;
workspace;
fontSize = 16;
% Multi-Scale Retinex on L Component of CIE-LAB Space
%I=imread("D:\Search Works\Underwater Dataset\raw-890 images\23_img_.png");
I=imread("23_img_.png");
Ilab = rgb2lab(I);
Ir = Ilab(:,:,1);
Ig = Ilab(:,:,2);
Ib = Ilab(:,:,3);
%%
%%%%%%%%%% Set the required parameters %%%%%%
G = 192;
b = -30;
alpha = 125;
beta = 46;
Ir_double=double(Ir);
%%
%%%%%%%%%% Set the Gaussian parameter %%%%%%
sigma_1=15; % Three Gaussian Kernels
sigma_2=80;
sigma_3=250;
[x, y]=meshgrid((-(size(Ir,2)-1)/2):(size(Ir,2)/2),(-(size(Ir,1)-1)/2):(size(Ir,1)/2));
gauss_1=exp(-(x.^2+y.^2)/(2*sigma_1*sigma_1)); %Calculate the Gaussian function
Gauss_1=gauss_1/sum(gauss_1(:)); %Normalization
gauss_2=exp(-(x.^2+y.^2)/(2*sigma_2*sigma_2));
Gauss_2=gauss_2/sum(gauss_2(:));
gauss_3=exp(-(x.^2+y.^2)/(2*sigma_3*sigma_3));
Gauss_3=gauss_3/sum(gauss_3(:));
%%
%%%%%%%%%%Operates on R component%%%%%%%
% MSR Section
Ir_log=log(Ir_double+1); % Converts an image to a logarithm field
f_Ir=fft2(Ir_double); % The image is Fourier transformed and converted to the frequency domain
%sigma = 15 processing results
fgauss=fft2(Gauss_1,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss); % Move the center of the frequency domain to zero
Rr=ifft2(fgauss.*f_Ir); % After convoluting, transform back into the airspace
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr1=Ir_log-Rr_log;
%sigma=80
fgauss=fft2(Gauss_2,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr2=Ir_log-Rr_log;
%sigma=250
fgauss=fft2(Gauss_3,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr3=Ir_log-Rr_log;
Rr=0.33*Rr1+0.34*Rr2+0.33*Rr3; % Weighted summation
MSR1 = Rr;
SSR1 = Rr2;
%SSR
min1 = min(min(SSR1))
max1 = max(max(SSR1))
% Rescale the SSR image into the range of L, which is 0-100.
SSR1 = rescale(SSR1, 0, 100);
%MSR
min2 = min(MSR1(:))
max2 = max(MSR1(:))
% Rescale the MSR image into the range of L, which is 0-100.
MSR1 = rescale(MSR1, 0, 100);
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Assign SSR1 to L and take Ig as A and Ib as B and create an LAB image.
ssr = cat(3,SSR1,Ig,Ib);
ssr = double(ssr);
ssr = lab2rgb(ssr); % Convert from lab back to RGB.
% Assign MSR1 to L and take Ig as A and Ib as B and create an LAB image.
msr = cat(3,MSR1,Ig,Ib);
msr = double(msr);
msr = lab2rgb(msr); % Convert from lab back to RGB.
subplot(3,2,1);
imshow(I, []);
impixelinfo;
title('Original') % Show the original image
subplot(3,2,2);
imshow(ssr, []);
impixelinfo;
title('SSR')
subplot(3,2,3);
imshow(msr, []);
impixelinfo;
title('MSR')
% Plot histograms of msr and ssr
subplot(3,2,4);
histogram(ssr, 500);
grid on;
title('Histogram of ssr (L) channel')
subplot(3,2,5);
histogram(msr, 500);
grid on;
title('Histogram of msr (L) channel')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021695/image.png)
0 Comments
More Answers (2)
DGM
on 4 Jun 2022
Edited: DGM
on 4 Jun 2022
The output of rgb2lab() is already of class double(). L is in the range of 0-100. Rescaling L to [0 255] and casting to uint8 will just result in nonsense output. I'm assuming the result was intended to be scaled to the full range of the L channel. If so, scale it to [0 100]. Bear in mind that significant changes to L in an image with significantly saturated colors will result in points moving out of gamut and potentially being truncated by subsequent operations, which may result in some color distortion. Whether that's important depends on the intent.
I = imread('peppers.png');
Ilab = rgb2lab(I);
Ir = Ilab(:,:,1); % these aren't R,G,B, so why call them that?
Ig = Ilab(:,:,2);
Ib = Ilab(:,:,3);
%%%%%%%%%%Set the required parameters%%%%%%
G = 192;
b = -30;
alpha = 125;
beta = 46;
Ir_double=double(Ir); % it's already double
%%%%%%%%%%Set the Gaussian parameter%%%%%%
sigma_1=15; %Three Gaussian Kernels
sigma_2=80;
sigma_3=250;
[x, y]=meshgrid((-(size(Ir,2)-1)/2):(size(Ir,2)/2),(-(size(Ir,1)-1)/2):(size(Ir,1)/2));
gauss_1=exp(-(x.^2+y.^2)/(2*sigma_1*sigma_1)); %Calculate the Gaussian function
Gauss_1=gauss_1/sum(gauss_1(:)); %Normalization
gauss_2=exp(-(x.^2+y.^2)/(2*sigma_2*sigma_2));
Gauss_2=gauss_2/sum(gauss_2(:));
gauss_3=exp(-(x.^2+y.^2)/(2*sigma_3*sigma_3));
Gauss_3=gauss_3/sum(gauss_3(:));
%%%%%%%%%%Operates on R component%%%%%%%
% MSR Section
Ir_log=log(Ir_double+1); %Converts an image to a logarithm field
f_Ir=fft2(Ir_double); %The image is Fourier transformed and converted to the frequency domain
%sigma = 15 processing results
fgauss=fft2(Gauss_1,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss); %Move the center of the frequency domain to zero
Rr=ifft2(fgauss.*f_Ir); %After convoluting, transform back into the airspace
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr1=Ir_log-Rr_log;
%sigma=80
fgauss=fft2(Gauss_2,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr2=Ir_log-Rr_log;
%sigma=250
fgauss=fft2(Gauss_3,size(Ir,1),size(Ir,2));
fgauss=fftshift(fgauss);
Rr= ifft2(fgauss.*f_Ir);
min1=min(min(Rr));
Rr_log= log(Rr - min1+1);
Rr3=Ir_log-Rr_log;
Rr=0.33*Rr1+0.34*Rr2+0.33*Rr3; %Weighted summation
MSR1 = Rr;
SSR1 = Rr2;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SSR
SSR1 = mat2gray(SSR1)*100;
ssr = cat(3,SSR1,Ig,Ib);
ssr = lab2rgb(ssr);
%MSR
MSR1 = mat2gray(MSR1)*100;
msr = cat(3,MSR1,Ig,Ib);
msr = lab2rgb(msr);
subplot(3,2,1);imshow(I);title('Original') %Show the original image
subplot(3,2,2);imshow(ssr);title('SSR')
subplot(3,2,3);imshow(msr);title('MSR')
% roughly 9% of pixels are OOG
nnz(ssr>1 | ssr<0)/numel(ssr)
nnz(msr>1 | msr<0)/numel(msr)
3 Comments
DGM
on 4 Jun 2022
Edited: DGM
on 4 Jun 2022
Let me preface this by saying that I made no attempt to troubleshoot your retinex implementation. If the output is grossly unexpected, then the issue probably is in the retinex code and mostly not due to truncation. That's just a guess.
It depends if truncation is actually a problem for your needs.
I give an explanation of what's happening geometrically in the answer here. That example discusses YCbCr, but the effect is the same in LAB. The result is that truncation due to a change in L will end up producing an unexpectedly smaller change in L and an unintended change to saturation and hue. Sometimes the distortion is unobjectionable visually, but if you're expecting a certain change in L, it might matter.
Consider the examples:
inpict = imread('peppers.png');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021395/image.png)
inlab = rgb2lab(inpict);
inlab(:,:,1) = inlab(:,:,1)+10;
outpict = lab2rgb(inlab); % about 5% oog
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021400/image.png)
The MIMT LCH tools can optionally truncate in LCH prior to conversion. This moves points radially, so L and H don't get distorted. OOG points are represented at the maximal chroma available for their given L and H prior to conversion.
% using MIMT lch tools
inlab = rgb2lch(inpict,'lab');
inlab(:,:,1) = inlab(:,:,1)+10;
outpict2 = lch2rgb(inlab,'lab','truncatelch');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1021410/image.png)
I don't know if it's possible to wrangle makecform() to do any sort of similar gamut mapping. I don't really ever use it.
See Also
Categories
Find more on Color in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!