Examples of Poisson noise and extraction of information 

Matlab codes 

% Program to improve S/N ratio in small factor steps of 2^0.1 = 1.0718 of a

% 16-bit tif photograph (copyright Nyah Willmott)

 

clear; close all;

vid = VideoWriter('poissonNoise.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 15;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

I = imread('zurich.tif'); % input image

J = rgb2gray(I); % convert to b&w

K = imresize(J,0.35); % rescale image to more practical size

for i = 16:-0.1:0

    factor = 2^i;

    L = K/factor; % reduce bit range. Reduces each pixel to round(val/factor)

    M = imnoise(L,'poisson'); % Add Poisson noise

    N = M*factor; % Multiply up again, now with noise added

    imshow(N);

    bits = round(16.4 - i);

    str1 = num2str(bits,'% 2i');

    if (bits >=1 )

        if (bits == 1)

            str2 = ' bit';

        else

            str2 = ' bits';

        end

        strTot = [str1,str2];

        hText = text(1600, 1190, strTot, 'FontSize',16,'Color','white');

    end

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);

% Program to improve S/N ratio in small factor steps of 2^0.1 = 1.0718 of a

% 16-bit tif photograph (copyright Nyah Willmott)

 

clear; close all;

vid = VideoWriter('poissonNoiseCroppedRegion.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 15;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

I = imread('zurich.tif','PixelRegion',{[1554,1680],[1760,1873]}); % input image

J = rgb2gray(I); % convert to b&w

for i = 16:-0.1:0

    factor = 2^i;

    L = J/factor; % reduce bit range. Reduces each pixel to round(val/factor)

    M = imnoise(L,'poisson');

    N = M*factor;

    imshow(N);

    bits = round(16.4 - i);

    str1 = num2str(bits,'% 2i');

    if (bits >=1 )

        if (bits == 1)

            str2 = ' bit';

        else

            str2 = ' bits';

        end

        strTot = [str1,str2];

        hText = text(1600, 1190, strTot, 'FontSize',16,'Color','white');

    end

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);

% Program to generate signal with noise and background with noise, fit them

% via a least squares fitting function provided by matlab, and plot the

% progress of the fit with number of scans.

% The function is a constant background with amplitude = 25, a peak signal

% with height = 2, centred at x = 250 and with a SD value of 16.

 

clear;

close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

vid = VideoWriter('snr.mp4','MPEG-4');

vid.FrameRate = 60;    % Default 30

vid.Quality = 100;    % Default 75

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

 

philBlue = [0.3 0.28 0.55];

 

numScans = 2000; % Number of scans

scanLength = 500; % Length of any given scan

sigPos = scanLength/2; % Place peak signal ideally at centre of scan

sigma = 16; % SD in Gaussian description of peak signal

xx = 1:scanLength; % x-coordinate of scan

bckav = 25; % Average background signal

sigpk = 2; % Average signal peak above background

backSig = zeros(1,scanLength); % Initialized background signal

sig = zeros(1,scanLength); % Initialized peak signal

 

% Initialize plots of progress in fitting

sigPkArray = zeros(1,numScans);

sigPosArray = zeros(1,numScans);

sigSigmaArray = zeros(1,numScans);

bckArray = zeros(1,numScans);

 

for i = 1:numScans

    subplot(4,5,[1,2,3,4,6,7,8,9,11,12,13,14,16,17,18,19])

    newplot

    hold off

    % Simulated background with Poisson noise

    backSig = backSig+poissrnd(bckav*ones(1,scanLength));

    % Idealized peak signal

    perfSig = sigpk*exp(-(xx - scanLength/2).^2/(2*sigma^2));

    % Simulated peak signal with Poisson noise

    sig = sig+poissrnd(perfSig);

    % Add up simulated noisy signals and normalized by dividing by i

    totSig = (sig + backSig)/i;

    plot(xx,totSig,'color',philBlue,'LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',18);

    ymin = 0.98*min(totSig);

    ymax = 1.02*max(totSig);

    ylim([ymin ymax]);

    fun = @(x,xx)x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);

    x0 = [10,scanLength/2.2,sigma*2,10]; % Random starting guess

    x = lsqcurvefit(fun,x0,xx,totSig);

    hold on

    fity = x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);

    plot(xx,fity,'color','red','LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',14);

    str1 = 'Scan ';

    str2 = num2str(i,'% 4.0f');

    strTot = [str1,str2];

    hTextYpos = ymin+(ymax-ymin)*0.88;

    hText = text(0.88*scanLength, hTextYpos, strTot, 'FontSize',18);

    sigPkArray(i) = x(1);

    sigPosArray(i) = x(2);

    sigSigmaArray(i) = x(3);

    bckArray(i) = x(4);

 

 

    %hold on

    subplot(4,5,[5])

    plot(sigPkArray(1:i),'color','red','LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',14);

    xlim([1 numScans]);

    ylim([sigpk*0.5 sigpk*1.5]);

    hText = text(numScans*0.8, sigpk*1.25, 'a_0', 'FontSize',14);

 

    subplot(4,5,[10])

    plot(sigPosArray(1:i),'color','red','LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',14);

    xlim([1 numScans]);

    ylim([sigPos*0.98 sigPos*1.02]);

    hText = text(numScans*0.8, sigPos*1.01, 'a_1', 'FontSize',14);

 

 

    subplot(4,5,[15])

    plot(sigSigmaArray(1:i),'color','red','LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',14);

    xlim([1 numScans]);

    ylim([sigma*0.5 sigma*1.5]);

    hText = text(numScans*0.8, sigma*1.25, 'a_2', 'FontSize',14);

 

    subplot(4,5,[20])

    plot(bckArray(1:i),'color','red','LineWidth', 2);

    set(gca,'linewidth',2,'fontsize',14);

    xlim([1 numScans]);

    ylim([bckav*0.98 bckav*1.02]);

    hText = text(numScans*0.8, bckav*1.01, 'a_3', 'FontSize',14);

 

    % Store the frame

    frame = getframe(gcf);

    writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);