How do I get this code give me RGB results?

I wrote this code for image holograms. But I want this to give results in RGB. It is giving in black and white. How to change it? This is my image.

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

Answer 1

[Image Analyst]

Your loop over rgb extracts and operates on one color channel at a time, so of course it's gray scale. It does each color channel in turn, one at a time. Any particular color channel (R, G, or B) is monochrome (gray scale).

[Diptayan Dasgupta]

Yes but isnt there any other way to combine those 3 data into 1 color image?

[Image Analyst]

To combine matrices into an RGB, 3-channel image

rgbImage = cat(3, r, g, b);

[Diptayan Dasgupta]

I used this but it's showing error.

if rgb==1
    Ii=Red;
elseif rgb==2
    Green=Ii;
elseif rgb==3
    Blue=Ii;
end
IMAC=cat(3,Red,Blue,Green);
imshow(IMAC)
title('Color image')

[Image Analyst]

IMAC needs to either be a floating point, real image in the range 0-1, or a uint8 image in the range 0-255.

IMAC = rescale(real(IMAC), 0, 1));  % To make it double
IMAC = uint8(255 * rescale(IMAC, 0, 255));  % To make it uint8