How do I change this code for RGB images?

I wrote this code for black and white images. But when I change this for RGB images, it doesnt work. This is the source image.

clear all, close all;
Ii=imread('BMB2.jpg');% 256*256 pixels 8bit
figure; imshow(Ii);
title('Object pattern')
axis off
Ii=double(Ii);
for rgb=1:3
Ii=Ii(:,:,rgb);
PH=rand([256,256]);
Ii=Ii.*exp(2i*pi*PH); % add a random phase on the object
M=512;
I=zeros(512);
I(128:383,128:383)=Ii; % zero padding
z=15; %(cm, distance) PLEASE CHANGE THE DISTANCE TO 1, 5, 15, ETC.
w=6500*10^-8; %(cm, wavelength)
delta=0.005; % cm, pixel size 50um
r=1:M;
c=1:M;
[C, R]=meshgrid(c, r);
% Forward propagation (650nm)
p=exp(-2i*pi*z.*((1/w)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
A0=fftshift(ifft2(fftshift(I)));
Az=A0.*p;
E=fftshift(fft2(fftshift(Az))); % 1st order of the hologram
% Reconstruction (650nm)
p=exp(-2i*pi*(-z).*((1/w)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
A1=fftshift(ifft2(fftshift(E)));
Az1=A1.*p;
R1=fftshift(fft2(fftshift(Az1)));
R1=(abs(R1)).^2;
figure; imshow(R1/max(max(R1)));
title('Reconstructed image(650nm)')
axis off
% Reconstruction (450nm~650nm)
dw=50;
IMA=zeros(512,512);
for g=0:40;
w2=(20000-dw*g)*10^-8; % reconstruction wavelength
E2=E.*exp(2i*pi*sind(10)*(w-w2)/w/w2.*R*delta);
% phase mismatch due to the wavelength shift
p=exp(-2i*pi*(-z).*((1/w2)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
Az2=ifft2(fftshift(E2)).*(fftshift(p));
R2=fftshift(fft2(Az2));
R=(abs(R2)).^2; % summation of all wavelengths
IMA=IMA+R2;
end
IMA=IMA/max(max(IMA));
figure; imshow(IMA)
title('Reconstructed image(white light)')
axis off
end

Answer 1

[KSSV]

In this line

Ii=Ii(:,:,rgb);

You are overwrititng the m*n*3 image variable into a 1D matrix. Thus you lost your rgb original image and when you are trying to access 2 channel you will get error.

Rename the original image Ii into some other variable. Like:

clear all, close all;
Ii=imread('BMB2.jpg');% 256*256 pixels 8bit
figure; imshow(Ii);
title('Object pattern')
axis off
Ii0=double(Ii);
for rgb=1:3
Ii=Ii0(:,:,rgb);
PH=rand([256,256]);
Ii=Ii.*exp(2i*pi*PH); % add a random phase on the object
M=512;
I=zeros(512);
I(128:383,128:383)=Ii; % zero padding
z=15; %(cm, distance) PLEASE CHANGE THE DISTANCE TO 1, 5, 15, ETC.
w=6500*10^-8; %(cm, wavelength)
delta=0.005; % cm, pixel size 50um
r=1:M;
c=1:M;
[C, R]=meshgrid(c, r);
% Forward propagation (650nm)
p=exp(-2i*pi*z.*((1/w)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
A0=fftshift(ifft2(fftshift(I)));
Az=A0.*p;
E=fftshift(fft2(fftshift(Az))); % 1st order of the hologram
% Reconstruction (650nm)
p=exp(-2i*pi*(-z).*((1/w)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
A1=fftshift(ifft2(fftshift(E)));
Az1=A1.*p;
R1=fftshift(fft2(fftshift(Az1)));
R1=(abs(R1)).^2;
figure; imshow(R1/max(max(R1)));
title('Reconstructed image(650nm)')
axis off
% Reconstruction (450nm~650nm)
dw=50;
IMA=zeros(512,512);
for g=0:40;
w2=(20000-dw*g)*10^-8; % reconstruction wavelength
E2=E.*exp(2i*pi*sind(10)*(w-w2)/w/w2.*R*delta);
% phase mismatch due to the wavelength shift
p=exp(-2i*pi*(-z).*((1/w2)^2-(1/M/delta)^2.*(C-M/2-1).^2-(1/M/delta)^2.*(R-M/2-1).^2).^0.5);
Az2=ifft2(fftshift(E2)).*(fftshift(p));
R2=fftshift(fft2(Az2));
R=(abs(R2)).^2; % summation of all wavelengths
IMA=IMA+R2;
end
IMA=IMA/max(max(IMA));
figure; imshow(IMA)
title('Reconstructed image(white light)')
axis off
end