Search code examples
imagematlabimage-processingastronomyimagefilter

Image Stacking for noise reduction using MATLAB

I'm new to image processing and started using MATLAB for Astrophotography processing. I'm trying to process 10 corrupted images (same image but mixed with different noise) of the planet Saturn using MATLAB. I learned that by stacking the 10 images together leads to a noise reduced picture with high PSNR and tried the below coding to make it work.

But the output looks like an unclear saturated image with no noise reduction.

enter image description here

Can you please look at the code below and show me where I went wrong?

%% We are going to stack the 10 corrupted images and finally calculate the PSNR SSIM
clearvars;% Clear all the variables
close all;

load('planetdata.mat'); %to load the corrupted Image set (4-D uint8)
Clean = imread('Clean Image of Saturn.jpg');%Clean Image of Saturn.600x800x3 uint8
planet1(: , :, :)   = planetdata(1, :, :, :);%One corrupted Image as reff

% Set the number of images to stack is 10
stack_number = 10;

% Lets use Clean image as reference of dimensions required
im_x = size(Clean, 1);
im_y = size(Clean, 2);
im_z = size(Clean, 3);

% Lets Generate a blank image for image stacking
resultIM = uint8(zeros(im_x, im_y, im_z));

% Iterate through the images to stack
for i = 1:1:stack_number

% Read in the target object image
 CorruptIM(: , :, :)   = planetdata(i, :, :, :);

% Perform image stacking using the target object image
 resultIM = resultIM + CorruptIM;

end

% resultIM = resultIM / stack_number; 

%% Lets Display Results
workspace;  % to Make sure the work space panel is showing.
fontSize = 15;
figure;
subplot(1, 3, 1);
imshow(Clean);
title('Clean Image', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Position', get(0,'Screensize')); 
% Give a name to the title bar.
set(gcf,'name','Stacking','numbertitle','off')

% Display one corrupt image as reference
subplot(1, 3, 2);
imshow(planet1);
title('Corrupt Image 1 : Ref', 'FontSize', fontSize);

% Display Stacked image
subplot(1, 3, 3);
imshow(resultIM);
title('Stacked Image', 'FontSize', fontSize);

%% PSNR AND SSIM Calculation
%Lets Find PSNR for For Resultant Image

[row,col]   = size(Clean);
size_host   = row*col;
o_double    = double(Clean);
w_double    = double(resultIM);

s=0;
for j = 1:size_host % the size of the original image

s = s+(w_double(j) - o_double(j))^2 ; 
end

mes     =s/size_host;
psnr    =10*log10((255)^2/mes);
fprintf('The PSNR value for Stacked Image is %0.4f.\n',psnr);


%Lets Find SSIM for resultant Image
[ssimval, ssimmap] = ssim(uint8(resultIM),Clean);
fprintf('The SSIM value for Stacked Image is %0.4f.\n',ssimval);
Solution

  • I think this one statement pretty much says it all (emphasis mine):

    But the output looks like an unclear saturated image with no noise reduction.

    It looks like your images are in fact saturating at the upper limit for a uint8 variable, which is the data type for your result image resultIM and your data matrix planetdata. As you keep adding images, the pixel values saturate at the maximum value of 255 for unsigned 8-bit integer types.

    What you'll need to do is convert to a larger data type for your intermediate calculations, such as a larger integer type or a floating point type. Then, once your calculations are complete (e.g. dividing the sum by the image stack size to get an average image), you can scale and/or round your data as needed and convert it back to a uint8 type.