Problem in getting proper output image format

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')

Answer 1

[Image Analyst]

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')

Answer 2

[Image Analyst]

It's probably floating point. So use [] in imshow

imshow(msr, [])

[D.Regan]

Thank You for your kind response. I tried in that way, got white blank image output only. Any other ways to solve?

Answer 3

[DGM]

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)

ans = 0.0888

nnz(msr>1 | msr<0)/numel(msr)

ans = 0.0886

[D.Regan]

Thank You. Then How can this out of gamut problem be solved?

[DGM]

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');

inlab = rgb2lab(inpict);
inlab(:,:,1) = inlab(:,:,1)+10;
outpict = lab2rgb(inlab); % about 5% oog

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');

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.

[D.Regan]

Thank You