0

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 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 a unclear saturated image with no noise reduction. Output Image

Can you please look at the code and point me where i went wrong.

Thank You

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

closed as off-topic by Olivier, mattdm, inkista, Philip Kendall, jwenting May 24 '17 at 7:06

  • This question does not appear to be about photography within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • 9
    I'm voting to close this question as off-topic because it is about programming, not photography. This question belongs on Stack Overflow. – scottbb May 8 '17 at 2:47
  • 3
    resultIM is an unsigned 8 bit array. You can't accumulate the stack images which are also unsigned 8 bits w/o overflowing and wrapping. – doug May 9 '17 at 2:23
  • 4
    While the programming part is likely off-topic, the problem is obvious from a photographic standpoint. I think this is answerable at a conceptual level if not a "fix my code" level. – AJ Henderson May 23 '17 at 13:43
  • @AJHenderson that's a good point. I would reverse my close vote, but that's not possible. Well, I reverse in principle. =) – scottbb May 23 '17 at 14:29
3

Image stacking is a process by which you can reduce noise, but it doesn't work by adding the images together additively, but rather averaging them. The reason that stacking works is that signal from the same photo taken multiple times will be the same, but random noise will be different each time.

If you average the images, the noise will tend to be reduces as it is not the same between images.

  • If you look at the code, it is an average of multiple images. – Olivier May 23 '17 at 17:14
  • @Olivier - I can't speak to the contents of the matlab code, but I can speak to the end result image. The end result image appears to be much brighter (probably either an additive or multiplicative blending). Therefore it can't be properly doing an average blend. How to specifically do an average image blend in matlab is really not on topic for the site. – AJ Henderson May 23 '17 at 17:18
  • 1
    You are right about something weird in the Matlab code. It seems to come from bad implementation : adding uint8 values multiple times into an uint8 one => instant overlow – Olivier May 23 '17 at 17:22
2

Your noise may not be zero mean, but the larger issue is that you are adding to a uint8 array. These values go up to 255, and will overflow likely by the second image.

You should do your intermediate steps in floating point, and then later cast to uint8 and perhaps normalize if necessary.

1

Yes, this answer should be in SO, along with the question. Here goes:

Your "noise" is almost certainly not zero-mean. You are adding all the noise values in all frames, which naturally leads to an average value for each of R,G,B that is the same, i.e. a white (or grey) pixel. Your final image is unhelpfully autoscaled by imshow . Since you failed to provide the code you used to synthesize noise, I might be wrong, but please review your noise-only files and try again.

edit: adding Sean's info:

A more immediate issue, however, is that because resultIM is of type uint8 it clips values over 255. By adding 10 different copies of CorruptIM, it looks like a large number of pixels are being turned white (255)

  • 1
    Excellent observation about non-zero mean noise. A more immediate issue, however, is that because resultIM is of type uint8 it clips values over 255. By adding 10 different copies of CorruptIM, it looks like a large number of pixels are being turned white (255). – Sean May 8 '17 at 21:35
  • 1
    @Sean thanks - I'd missed that class assignment. As a side note, it drives me batty that MATLAB coerces all objects in an expression to the lowest class ( uint8 + double is coerced to uint8). :-( – Carl Witthoft May 9 '17 at 11:24

Not the answer you're looking for? Browse other questions tagged or ask your own question.