Recolouring Sentinel2 multispectral image

Hello,

I am trying to recolour an image composed of the red, green and blue channels of Sentinel 2 imagery (Bands 4, 3 and 2, respectively). I have processed this imagery in the SNAP toolbox and exported in netCDF4-BEAM format.

Firstly, I read in the required image bands and create the RGB image:

address = %file location on computer;
ncid = netcdf.open(address);

bands{1} = transpose(im2double(netcdf.getVar(ncid,4)));    %Red channel
bands{2} = transpose(im2double(netcdf.getVar(ncid,3)));    %Green channel
bands{3} = transpose(im2double(netcdf.getVar(ncid,2)));    %Blue channel

Image = cat(3,bands{1},bands{2},bands{3});    %RGB image, datatype of double

When I try to visualise the image with imshow(Image), the result is near monochromatic. I understand this is due to correlation between channels and so I wish to recolour the image.

What is the best way to do this?

I have tried to use the code below, however this makes no difference:

stretchedImage = imadjust(Image,stretchlim(Image));

I have also investigated the image histograms however I don't understand why there is a difference between histogram(Image) and imhist(Image). I think the differences in the outputs of these may be related to my issues with strethlim, hopefully you can help shine a light on this for me!

Many thanks.

[Subhadeep Koley]

Can you share your netCDF4-BEAM file? or any other example netCDF4-BEAM file, which has the same dimension as yours.

[Aidan Wood]

Here is a link to one of the files I have been using: https://drive.google.com/open?id=1p-EM1ledycSIJo2TaQJdDR0C-NNGorVM

Answer 1

[Subhadeep Koley]

Hi, 3D scatter plot on your image produces a linear trend of the Red-Green-Blue channel, indicating a high correlation among them.

This is the reason why your image looks near monochromatic

You can equalize the individual band histograms before concatenation to view the image in RGB color. Also, rescaling the image in 0-255 range will help.

Refer the code below.

close all; clc;
address = 'S2A_MSIL1C_20181215T135101_N0207_R024_T22LCR_20181215T153027_s2resampled.nc';
ncid = netcdf.open(address);
bands{1} = transpose(im2double(netcdf.getVar(ncid,4)));    % Red channel
bands{2} = transpose(im2double(netcdf.getVar(ncid,3)));    % Green channel
bands{3} = transpose(im2double(netcdf.getVar(ncid,2)));    % Blue channel

bands{1} = imadjust(bands{1}, stretchlim(bands{1}));       % Adjustment on Red channel
bands{2} = imadjust(bands{2}, stretchlim(bands{2}));       % Adjustment on Green channel
bands{3} = imadjust(bands{3}, stretchlim(bands{3}));       % Adjustment on Blue channel

bands{1} = histeq(bands{1});                               % histeq on Red channel
bands{2} = histeq(bands{2});                               % histeq on Green channel
bands{3} = histeq(bands{3});                               % histeq on Blue channel

RGBImage = cat(3, bands{1}, bands{2}, bands{3});           % RGB image, datatype of double
RGBImage = uint8(rescale(RGBImage, 0, 255));               % Typecasted to uint8 and rescaled to 0-255 range for visualization
figure; imshow(RGBImage, []);

Otherwise, you can use the function decorrstretch also.

close all; clc;
address = 'S2A_MSIL1C_20181215T135101_N0207_R024_T22LCR_20181215T153027_s2resampled.nc';
ncid = netcdf.open(address);
bands{1} = transpose(im2double(netcdf.getVar(ncid,4)));    % Red channel
bands{2} = transpose(im2double(netcdf.getVar(ncid,3)));    % Green channel
bands{3} = transpose(im2double(netcdf.getVar(ncid,2)));    % Blue channel

RGBImage = cat(3, bands{1}, bands{2}, bands{3});           % RGB image, datatype of double
RGBImage = uint8(rescale(RGBImage, 0, 255));               % Typecasted to uint8 and rescaled to 0-255 range for visualization

RGBImage = decorrstretch(RGBImage, 'Tol', 0.03);           % decorrstretch on the RGBImage

figure; imshow(RGBImage, []);

3D scatter plot on your image after enhancement looks like below,

Hope this helps!

[Aidan Wood]

Thank you for your reply - it has been very helpful!

As part of my project I am standardising each image channel before applying PCA. When I try to visualise the standardised image (mean = 0, s.t.d = 1) I was not getting the same results, but by changing the second to final line to

RGBimage = im2uint8(RGBimage);

I get the same result.

If anyone is interested in replicating the RGB image respresentation in the SNAP toolbox, I have included the code I worked out below:

close all; clc;
address = 'S2A_MSIL1C_20181215T135101_N0207_R024_T22LCR_20181215T153027_s2resampled.nc';
ncid = netcdf.open(address);
R = transpose(im2double(netcdf.getVar(ncid,4)));    % Red channel
G = transpose(im2double(netcdf.getVar(ncid,3)));    % Green channel
B = transpose(im2double(netcdf.getVar(ncid,2)));    % Blue channel

R = rescale(R,0,1);
G = rescale(G,0,1);
B = rescale(B,0,1);

rIndices = find(abs(R)>0.2);
gIndices = find(abs(G)>0.2);
bIndices = find(abs(B)>0.2);

R(rIndices) = NaN;
G(gIndices) = NaN;
B(bIndices) = NaN;

RGB = cat(3, R, G, B);

RGB = uint8(rescale(RGB, 0, 255));

imshow(RGB, [])

[Aidan Wood]

Thanks also for the expanded answer!